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We revisit the model of a quantum Brownian oscillator linearly coupled to an environment of 
quantum oscillators at finite temperature. By introducing a compact and particularly well-suited 
formulation, we give a rather quick and direct derivation of the master equation and its solutions 
for general spectral functions and arbitrary temperatures. The flexibility of our approach allows for 
an immediate generalization to cases with an external force and with an arbitrary number of Brow- 
nian oscillators. More importantly, we point out an important mathematical subtlety concerning 
boundary-value problems for integro-differential equations which led to incorrect master equation 
coefficients and impacts on the description of nonlocal dissipation effects in all earlier derivations. 
Furthermore, we provide explicit, exact analytical results for the master equation coefficients and 
its solutions in a wide variety of cases, including ohmic, sub-ohmic and supra-ohmic environments 
with a finite cut-off. 
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I. INTRODUCTION 
A. New Results placed in Background Context 

An open quantum systems (OQS) [ ] refers to a quan- 
tum system interacting with an environment, which could 
be multi-partite, possessing many more degrees of free- 
dom (it could also be identified as the remaining "irrel- 
evant" degrees of freedom of the system itself). An en- 
vironment in some simplified modeling can be described 
in terms of its spectral density and parametrized by its 
temperature. Its influence on the open system can be 
expressed in terms of fluctuations (vacuum and thermal) 
and noises (the most general form can be colored and 
multiplicative). A theory of OQS describes the nature 
and dynamics of this system as a result of such interac- 
tions, which manifest in quantum dissipation and diffu- 
sion and can alter significantly the quantum coherence, 
entanglement and correlation properties of the otherwise 
closed quantum system. The familiar quantum statistical 
mechanics is the extreme limiting case when the system 
remains in equilibrium through interaction with a ther- 
mal or chemical reservoir. 

Open quantum system is the theoretical construct suit- 
able for the investigation of the properties and dynamics 
of nonequilibrium quantum systems in the Langevin vein 
(as distinguished from the Boltzmann vein, which consid- 
ers closed systems albeit often with a hierarchical struc- 
ture; see, e.g., Ref. [2]). It plays an important role in ad- 
dressing the fundamental issues such as the quantum-to- 
classical transition through the environment-induced de- 
coherence mechanism [3, 4]. For practical purposes it has 
been effectively applied to exciting phenomena in many 
new directions of micro- and meso-physics in the last 
two decades, made possible by innovative experiments 
aided by technological advances in high-precision instru- 
mentation. These include the areas of superconductiv- 
ity such as quantum dissipative tunneling in SQUIDs [ - 
], atomic and quantum optical systems using ultrafast 
lasers with atoms in cavities and optical lattices [8-10], 
as well as nanoelectromechanical devices [11, 12] which 
have great potential in physical, chemical and bioscience 
applications. For an accurate description of the system's 
properties and evolution in these processes, the effects of 
its interaction with the environment are essential. 



Quantum Brownian motion (QBM) of an oscillator 
coupled to a thermal bath of quantum oscillators has 
been extensively studied as a canonical model for open 
quantum systems because there is a considerable amount 
of insight that one can learn from it while being treat- 
able analytically to a significant degree. In this paper we 
continue the lineage of work on QBM via the influence 
functional path-integral method of Feynman and Vernon 
[13] used by Caldeira and Leggett [ ] to derive a mas- 
ter equation for a high-temperature ohmic environment, 
which corresponds to the Markovian regime. Following 
this, Caldeira, Cerdeira and Ramaswamy [ ] derived the 
Markovian master equation for the system with weak 
coupling to an ohmic bath, which was claimed to be valid 
at arbitrary temperature (see Sec. V C for a critique of 
this claim). At the same time Unruh and Zurek [16] de- 
rived a more complete and general master equation that 
incorporated a colored noise at finite temperature, but 
there is a problem with their fluctuation-dissipation re- 
lation (see Ref. [ ]). Finally, in a path-integral calcula- 
tion from first principles, Hu, Paz and Zhang (HPZ) [17] 
derived a master equation for a general environment (ar- 
bitrary temperature and spectral density), barring cer- 
tain subtle errors in the coefficients, which lead to in- 
accurate treatment of the nonlocal dissipation cases, as 
we will discuss. After that, this equation has been red- 
erived by a number of authors. Halliwell and Yu [18] ex- 
ploited the phase-space transformation properties of the 
Wigner function for the full system plus environment and 
derived a Fokker-Planck equation corresponding to the 
HPZ equation. Calzetta, Roura and Verdaguer (CRV) 
[19, 20] derived it using a stochastic description for open 
quantum systems based on Langevin equations, whereas 
Ford and O'Connell [21] employed a somewhat related 
method via the quantum Langevin equation [22] and ob- 
tained also the solution to the HPZ equation for a Gaus- 
sian wave-packet. 

The present paper's contribution to this legacy is three- 
fold: 

1. We have completely determined the precise form of 
the HPZ master equation coefficients and pointed 
out a problem with earlier derivations for nonlocal 
dissipation (Sec. IIIB). 

2. We have found concise and efficient solutions to the 
master equation with a number of exact nonpertru- 
bative analytical results (Sec. IV). 

3. We have extended the theory to that of a system 
of multiple oscillators bilinearly coupled amongst 
themselves and to the bath in an arbitrary fashion 
while acted upon by classical forces (Sec. VII). 

In this paper we will follow the approach of CRV in 
Refs. [19, 20] and make use of a stochastic description 
whose central element is a Langevin equation for the 
dynamics of the open quantum system. This offers an 
efficient mathematical tool for obtaining all the quan- 
tum properties of the system. An important feature of 
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the present approach is the reformulation in phase-space 
(rather than position space) together with the use of vec- 
tor and matrix notation. The combination of all these 
elements makes this new approach far more flexible and 
compact. For example, we are able to derive the general 
expression for the solution of the master equation in es- 
sentially two short lines [see Eq. (36)]. The flexibility of 
our formalism is also illustrated by the straightforward 
generalizations to the cases of an external force (this is 
nontrivial for nonlocal dissipation) and an arbitrary num- 
ber of system oscillators that will be presented. This 
goes far beyond previous generalizations of the theory 
[ !] which assume specific forms of coupling. 

One of our key contributions, however, is uncovering 
a significant shortcoming of earlier results for the master 
equation coefficients. We point out a subtlety involving 
boundary conditions for solutions of integro-differential 
equations and explain how certain properties that hold 
for ordinary differential equations are not true for non- 
local dissipation. These properties had always been em- 
ployed erroneously, in one way or another, when deriving 
the expressions for the master equation coefficients, even 
those which were then evaluated numerically. This long- 
standing error could have deep implications for regimes 
where the effects of nonlocal dissipation are significant 
and one should be cautious with all results for those cases 
reported in the literature. 

Taking into account the aspect mentioned in the pre- 
vious paragraph, and using our compact formulation, we 
have provided a relatively simplified expression for the 
correct master equation. Moreover, one can also ob- 
tain the general solution to the master equation in terms 
of the matrix propagator of a linear integro-differential 
equation, and see that at late times it tends to a Gaussian 
state completely characterized by a constant covariancc 
matrix. For odd meromorphic spectral functions, and 
many others, we are able to reduce the calculation of this 
covariance matrix to a simple contour integral and obtain 
exact nonperturbative results for finite cut-off and arbi- 
trarily strong coupling. This includes examples of ohmic, 
sub-ohmic and supra-ohmic environments; and from this 
late-time covariance one can immediately obtain the late- 
time diffusion coefficients as well. Our results generalize 
the work of Anastopoulos and Halliwell [24] as well as 
Ford and O'Connell [21], who already found the late time 
state to be a Gaussian, and the earlier work of Hu and 
Zhang [25, 26] on the generalized uncertainty function for 
Gaussian states. 

In addition, working with Laplace transforms and then 
transforming back to time domain, we manage to find 
the exact solutions for the propagators associated with 
the integro-differential equations corresponding to ohmic, 
sub-ohmic and supra-ohmic environments with a finite 
cut-off. This enables us to gain very valuable informa- 
tion on the dynamics of the system. For instance, for 
an ohmic environment one can show that using the local 
approximation for the propagator is a valid approxima- 
tion in the large cut-off limit, which makes it possible 



to obtain relatively manageable analytic results for the 
diffusion coefficients at all times. Furthermore, the exact 
solution of a specific sub-ohmic environment reveals that 
long-time correlations (due to excessive coupling with IR 
modes of the environment) give rise to contributions to 
the propagator that decay at late times like power laws. 
This invalidates the use of an effectively local descrip- 
tion at late times, whose contributions decay exponen- 
tially, and provides a clear example of a situation where 
nonlocal dissipation needs to be properly dealt with. Fi- 
nally, studying the exact solutions for some particular 
supra-ohmic environment we also find significant nonlo- 
cal effects which are due in this case to the UV regulator 
function. This leads to a marked cut-off sensitivity of the 
momentum covariance that had not been noticed before. 



B. Key Points and Organization 

Those readers who want to find out quickly the prob- 
lem with earlier derivations of the master equation can 
simply read Sec. II to get acquainted with our nota- 
tion and formalism and go to Sec. IIIB, where the mas- 
ter equation is derived, aided perhaps by D, which ex- 
plains in detail the key mathematical subtlety concerning 
integro-differential equations and its implications for the 
existing derivations. They may also find Sec. VI valuable 
since it contains specific examples where nonlocal dissi- 
pation effects give dominant contributions and can lead 
to significant discrepancies from previous results. 

The other useful results are mentioned below alongside 
a description of how this paper is organized. 

The key framework providing the stochastic descrip- 
tion for an open quantum system in terms of a Langevin 
equation and its compact phase-space formulation is in- 
troduced in Sec. II, where a very simple derivation of 
the general solution for the state evolution of the sys- 
tem is given. The problems with previous derivations 
are pointed out and the correct derivation of the master 
equation is given in Sec. III. The master equation is then 
solved using the method of characteristic curves and the 
solution is shown to be equivalent to that obtained in a 
more straightforward manner from the Langevin equa- 
tion. 

The general solution of the master equation is em- 
ployed in Sec. IV to discuss general properties of the state 
evolution of the QBM subsystem, tending to a Gaussian 
stationary state at late times. A very simple and intuitive 
picture of environment-induced decoherence in terms of 
the reduced Wigner function can be directly extracted, 
which could easily be made quantitative and precise. In 
addition, a generic discussion of late-time dynamics is 
provided. 

In Sec. V we find the exact nonlocal propagator for an 
ohmic environment with finite cut-off and identify a new 
regime at ultra-strong coupling. We provide exact non- 
pertrubative results for the late-time thermal covariance 
and full-time results for the diffusion coefficients in the 
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large cut-off limit. 

Explicit examples of sub-ohmic and supra-ohmic spec- 
tral functions arc considered in Sec. VI for which the ex- 
act propagator is computed and dominant contributions 
from nonlocal dissipation effects are found (of IR origin 
in one case and UV in the other). 

The generalization to a system of multiple oscillators 
bilinearly coupled to themselves and the bath in arbitrary 
fashion and acted upon by classical forces is presented in 
Sec. VII. Finally, in Sec. VIII we summarize our results 
and discuss their main implications as well as possible 
applications. 

In addition to a couple of appendices on special func- 
tions and properties of Laplace transforms for reference 
purposes, C contains technical aspects concerning diver- 
gences of the dissipation kernel and frequency renormal- 
ization, as well as initial kicks and a discussion of diver- 
gences associated with uncorrelated initial states. 

D contains a detailed explanation of the mathematical 
subtlety involving boundary-value problems for integro- 
differential equations and a discussion of how it affected 
different classes of earlier derivations of the master equa- 
tions. The important formula for the late-time covariance 
in terms of a single frequency integral is derived in E, and 
the explicit analytic results for the diffusion coefficients 
of an ohmic environment at all times in the large cut-off 
limit are computed in F. 

Throughout the paper we use units with h = Jcb = !• 



II. THE LANGEVIN EQUATION 

A. General Theory 

The Lagrangian of a closed system consisting of a 
quantum Brownian oscillator with mass M, natural fre- 
quency f2 and coordinate x, bilinearly coupled with cou- 
pling constants c„ to an environment consisting of oscilla- 
tors with mass m n , natural frequency Lo n and coordinates 
x n , is most straightforwardly given by 



1 



M{x 2 



"bare"'' 



(x n 0J n x n j 



(1) 



E 



c n xx r . 



One introduces a "bare" frequency fibare because the 
interaction with the environment shifts the coefficient 
of the potential term by a certain amount <5f2 2 , given 
by Eq. (C3), so that the square of the actual fre- 
quency characterizing the subsystem of interest is given 



by Sli 



SCI . Alternatively, one can consider the fol- 



lowing Lagrangian: 



-M (x 2 
2 V 



2 J2\ 



Wx 



(2) 



E 



c„O s (t) 



where f2 corresponds to the actual frequency of the Brow- 
nian oscillator. For s (t) — 1 and provided that one 
identifies £1 2 with O^nre — * ms new Lagrangian is 
equivalent to that of Eq. (I) (further details on frequency 
renormalization and related issues are provided in C). In 
addition, we included a switch-on function 9 s (t) which 
vanishes at the initial time and smoothly increases to 
reach a constant unit value after a characteristic time- 
scale t s . While we consider initially uncorrelated states 
for the Brownian oscillator and the environment through- 
out the paper, which can sometimes lead to certain un- 
physical results, introducing a smooth switch-on function 
provides a way of effectively generating well-behaved ini- 
tial states with the high-frequency modes of the environ- 
ment properly correlated with the Brownian oscillator. 
Further discussion on this point can be found in C 2, but 
throughout the rest of the paper we will take 6> s (i) = I 
(or, equivalently, t s — 0) unless stated otherwise, and will 
only occasionally describe how the results would differ for 
a non-vanishing switch-on time. 



The subsystem corresponding to the quantum Brown- 
ian oscillator constitutes an open quantum system: while 
the evolution of the whole closed system is unitary, the 
Brownian oscillator (referred to as the "system" from 
now on) evolves non-unitarily due to the entanglement 
generated by the interaction with the environment. An 
important object characterizing the open system is the 
reduced density matrix, which results from taking the 
density matrix of the closed system and tracing out the 
environment: p T — Tr^p. The expectation value of ob- 
servables O that only depend on the system variables 
and are local in time can be directly obtained from it: 
(0)(t) = Tr [O pi(t)]. Given the density matrix for a 
continuous degree of freedom in position representation, 
one can always define the corresponding Wigner function: 



W T (X,p,t) 



1 f +oa / A 



A 



(3) 

which contains the same amount of information. See for 
instance Ref. [27] for a detailed description of the main 
properties of Wigner functions. In addition, the so-called 
dissipation and noise kernels (which involve respectively 
the commutator and anticommutator of the environment 
position operators in interaction picture) play an im- 
portant role when studying the open system dynamics 
[28, 29]. The case of a time-dependent coupling has been 
considered by Hu and Matacz [30], wherein all param- 
eters of the system and bath oscillators and their cou- 
plings were allowed to be time-dependent. When only 
the system-environment coupling is time-dependent, as 
in our case, and the initial state of the environment is 
a thermal state with temperature T, the dissipation and 
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noise kernels are given respectively by 

H(t,r) = - / dwsin[w(t-r)] I(u)0 e (t)0 a (r), (4) 



v(t, T ) = +j d " coth (^) cos[«(t-T)] J(w) e s {t)e s {r), 

(5) 



where /(w) is the spectral density function defined by 



^) = E^(— n). 



(6) 



It is often taken to be ohmic, i.e. = (2/ir)MjoU), 

but with a cut-off regulator so that it vanishes (or de- 
cays sufficiently fast) above some high-frequency scale A. 
However, more general spectral functions have been con- 
sidered before and will be considered here as well. 

It was shown in Ref. [19] that the quantum proper- 
ties of this kind of open systems can be entirely studied 
using a stochastic description whose central element is a 
Langevin equation of the form (L-x)(t) — £(i), where £(t) 
is a Gaussian stochastic source with a vanishing mean 
and correlation function equal to the noise kernel, i.e. 

= and (£(i)£( T )) ? = v{t,r). The dissipation 
kernel in turn appears in the Langevin integro-differential 
operator L, which is defined by 



(L ■ x)(t) = Mx{t) + Mn 2 x{t) + 2 dr /i(t, r) x(t) 

Jo 

+ M6n 2 9 2 (t)x(t), (7) 

and where <5f2 2 is given by Eq. (C3). One can then ex- 
press the time-evolving reduced Wigner function in terms 
of solutions of the Langevin equation and a double av- 
erage over their initial conditions, weighed with the re- 
duced Wigner function at the initial time, and over the 
realizations of the stochastic source [see Eq. (29) below]. 
Furthermore, one can also obtain the quantum correla- 
tion functions for system observables at multiple times 
(which in general cannot be obtained from the reduced 
Wigner function and its evolution via the master equa- 
tion) in terms of the solutions of the Langevin equation 
[19], as briefly illustrated in Sec. II D. See also Ref. [22] 
for a similar formulation involving a Langevin equation 
for operators in the Heisenberg picture. 

If we take a vanishing switch-on time, which amounts 
to discarding the switch-on function entirely, both the 
noise and dissipation kernels become time-translation in- 
variant. Moreover, it is convenient to introduce a damp- 
ing kernel "i{t—r) which is related to the dissipation kernel 
by fi(t, t) — (Jt(t—r) — M(d/dt)j(t—r) and is hence given 
by 

1 f°° I(u>) 
7 (i, r) = 7 (t-r) = — dw cob[w(*-t)] (8) 

M J LJ 

Note that this kernel is symmetric and positive definite 
like the noise kernel. Integrating by parts, the left-hand 



side of the Langevin equation can be written as follows 
(see C for further details): 

(L ■ x)(t) = Mx(t) + 2M [ drj(t-T) x(t) + Mfi 2 z(i) 
Jo 

+ 2M7(i)a(0), (9) 

The damping-kernel representation provides a cancela- 
tion of the frequency renormalization while introducing 
a slip in the initial conditions. This is caused by the 
last term on the right-hand side of Eq. (9), which cor- 
responds to a transient driving term proportional to the 
position of the system at the initial time. Leaving the 
slip term aside, one can show that all the (accumulated) 
energy dissipated through the nonlocal damping kernel 
term will be strictly positive (no amplification) as a con- 
sequence of the damping kernel being positive-definite. 



B. Solutions of the Langevin Equation 

The Langevin equation can be written as 

L ■ x = M (x + 27* x + tt 2 x) + 2Mx a j = f, (10) 

where * denotes the Laplace convolution, i.e. (A*B)(t) = 
J dTA(t—r)B{r), and Xq is the initial condition at t = 0. 
It is, thus, convenient to perform a Laplace transform 

/•OO 

f(s)=£{f}(s)= / dte- st f(t), (11) 
Jo 

under which the equation becomes purely algebraic. The 
Laplace transform of Eq. (10) is given by 

M (s 2 + 2s 7 (s) + n 2 ) £{s) = M (sx + x Q ) + i(s), (12) 
whose solution is 

x(s) = M (sx Q + ±0) G(s) + G(s)i(s), (13) 
1/M 



G(s) = 



(14) 



where terms proportional to the initial conditions xq and 
±0 correspond to the homogeneous solution while the 
noise term corresponds to the driven solution. G(t) satis- 
fies the initial boundary conditions G(0) = 0, G(0) = jL 
and fully determines the retarded Green function or prop- 
agator. In the time domain, the solution can be expressed 
as 

x{t) = M (x G(t) + x G(t)) + (G * £)(*)■ (15) 



1. Meromorphic Spectra 

For an ohmic environment in the infinite cut-off limit 
one has 7 (s) = 7o . More realistically, 7(5) will decay 
sufficiently fast at high s, implying a certain degree of 
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nonlocal dissipation (non-polynomial behavior in Laplace 
space). Thus, as illustrated by this example, one will gen- 
erally need to deal with non-polynomial damping kernels 
7(s). If 7(s) is a meromorphic function (i.e. analytic 
except for an isolated set of poles) , obtaining the inverse 
Laplace transform of G(s) amounts to calculating a sim- 
ple contour integral. 

On the other hand, given expression (8) for the damp- 
ing kernel, one can easily compute its Laplace transform: 



70) 



M J uj uj 2 + s 2 



(16) 



If we take the odd extension of the spectral density for 
negative frequencies, i.e. I(— \uj\) = — I(\lj\), then the 
integral can be recast as 



7(s) 



1 

2M 



UJ UJ 2 + s 2 



(17) 



which can be easily evaluated if the odd extension of I(uj) 
is meromorphic, e.g. for I(uj) ~ uj but not I(uj) ~ uj 2 . 
This is still less than ideal as the difficulty of solving the 
Langevin equation is more directly determined by the na- 
ture of the damping kernel. One would rather make the 
choice of damping kernel first (preferably in the Laplace 
domain) than derive it from the spectral density. Nev- 
ertheless, since the spectral density is still required to 
compute the noise kernel, we need the inverse relation- 
ship. Furthermore, as shown below, not every j(s) (even 
sufficiently regular ones) can be obtained from a spectral 
function through Eq. (17). 

Fortunately, Eq. (8) implies a simple relation between 
the spectral density and the Fourier transform of the 
damping kernel: I(uj) = — ^7(0;), and using Eq. (B14) 
applied to 7(0;) we get the following result for I(uj) in 
terms of the Laplace transform of the damping kernel: 



I(uj) = —Muj lira [^(e+tuj) + 7(e — uS)\ 

IT e->0 



(18) 



From this we see that meromorphic damping kernels re- 
sult in spectral densities which are odd meromorphic 
functions. Conversely, we have also seen that odd mero- 
morphic spectral densities lead to a meromorphic damp- 
ing kernel in Laplace space that can be obtained via con- 
tour integration through Eq. (17). We will thus refer 
to this class of odd meromorphic spectral densities and 
corresponding damping kernels as meromorphic spectra. 
Moreover, as we will see in later sections, given that 
Bromwich's formula for the inverse Laplace transform 
can also be computed as a contour integral, all the im- 
portant quantities for these meromorphic spectra are cal- 
culable via contour integration. 

Note that, as mentioned above, not every meromorphic 
function j(s) corresponds to a damping kernel that can 
be obtained from a spectral function through Eq. (17). 
This point can be seen by realizing that according to 
Eq. (18) different j(s) will give rise to the same spectral 
density as long as ~f(e+iuj)+j(e — iuj) is the same. Hence, 



if one wants to consider a candidate function 7(s), one 
should proceed as follows. Eq. (18) is first used to ob- 
tain the spectral density, which is then substituted into 
Eq. (17). If the initial candidate is recovered, it was a 
satisfactory one to begin with, otherwise it should be dis- 
carded, but the new damping kernel obtained in the last 
step is a valid one, which can be used instead. 



2. Phase-Space Representation 

If we introduce the phase-space coordinates z T = 
(x,p), the Langevin equation (10), together with the rela- 
tion p = mi, can be recast as a first-order linear integro- 
differential system of equations: 



H * z 



(19) 



where we introduced the boldface notation for vectors 
and matrices, £ T = (0, £) and the time-nonlocal pseudo- 
Hamiltonian H(t, t) = H(t— r) is given by 



H(r) 







y(r) 



M 



MQ 2 5(t) 2 7 (t) 



(20) 



Performing the Laplace transform of Eq. (19), which be- 
comes a purely algebraic equation, and rearranging the 
terms to express the solution in terms of the initial con- 
ditions and the stochastic source, one gets 



z(s) = *(s)z + *(s)£(s), 

MsG(s) G(s) 
M 2 s 2 G(s) - M MsG(s) 



*(*) 



(21) 
(22) 



where G(s) is the same propagator derived in the posi- 
tion representation and given by Eq. (14). Transforming 
back to the time domain, we can express the initial- value 
solutions as 



z(t) = *(t) z„ + (**£)(*). 
MG(t) G(t) 
M 2 G(t) MG(t) 



*(*) 



(23) 
(24) 



and <fr(i) can be identified as the matrix propagator as- 
sociated with the phase-space version of the Langevin 
equation, Eq. (19). 

Combining the result for z(t) as given by Eq. (23) with 
an analogous expression for the solution z(r) evaluated 
at an earlier time r < t, one can write z(r) in terms of 
z(t) and the stochastic source as follows: 

z(r) = *(r, t) z(t) - J dr 1 *(r, t) ^(t-r 1 ) £(r') 

- fdr' [*(r, t) *(i-r') - *(r-r')] Z(r'p5) 
Jo 

where we introduced the transition matrix <f?(t, r), which 
is defined as 



*(t,r) = *(i)* _1 (r) 



(26) 
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Note that 3?(t, r) ^ 4>(t — r) unless one has local dissi- 
pation. Thus, in the general case of nonlocal dissipation 
the last term on the right-hand side of Eq. (25) does 
not vanish and z(r) also depends on £(r') with t' < r. 
This means that, unlike with ordinary differential equa- 
tions, when boundary conditions z(i) are specified at a 
final time t, there is no truly advanced propagator for 
the inhomogeneous solutions of the integro-differential 
equation. One can still express the solution of such a 
final-value problem in terms of a matrix propagator (or 
Green's function in position space) with the right bound- 
ary conditions: 

z(r) = *(r, t) z(t) + fdr' # f (r, r') £(r'), (27) 
Jo 

where 

* f ( T , t') = -*(r, i) *(t-r') + 0(r-r') *(r-r'), (28) 

but one only has 3?f (r, r') = for r > t' in the case of 
strictly local dissipation. 

Such mathematical subtleties of final-value problems 
for integro-differential equations have been missed in the 
existing literature on the derivation of the master equa- 
tion for QBM models and could lead to significant dis- 
crepancies whenever the nonlocal effects of dissipation 
are important. A detailed discussion of this and related 
points is provided in D. 



C. Evolution of States 

As found in Ref. [19], the reduced Wigner function can 
be expressed in terms of the solutions of the Langevin 
equation and a double average over their initial condi- 
tions and the realizations of the stochastic source. Using 
the vector notation for phase-space variables introduced 
in the previous subsection, the result can be written as 



W I (z,t) = ((6(z(t)-z)). 



(29) 



with the averages over the initial conditions and the 
stochastic source defined as follows: 



^ J d*>--W T (x,0), 



(■■•>£ 



yj2ir det(z7j 



D£ • • • e 



(30) 

M-"- 1 -^ (31) 



where the right-hand side of Eq. (31) corresponds to the 
functional integral associated with the Gaussian stochas- 
tic source. The characteristic function of the Wigner 
function, regarded as a phase-space distribution, is given 
by its Fourier transform and it can be shown to take a 



rather simple form: 

W I (k,t)= [dze-^*((5{z-z(t)])J € , (32) 

e-^W) ) , (33) 
'*o/ e 

-ik T #(t)z \ / e -ik T (**£)(t)\ ^ (3^ 
/z \ '£ 

= W r (* T (i)k,0)e-5 kT -T(t)k ; (35) 
where the thermal covariance matrix <xr(i) is given by 

tr T (t)= I dr [ dr' 3>[t-T)v{T,T')$ T {t-T'), (36) 
Jo Jo 



u(t,t') 





v{t,t') 



(37) 



In the third equality above we used the initial-value so- 
lution (23) for z(i) to get Eq. (34), and in the last step 
we completed the square to calculate the Gaussian func- 
tional integral corresponding to the noise average in order 
to obtain the final result in Eq. (35). Note that for our 
Lagrangian, the stochastic force £ only has a momen- 
tum component and, therefore, all the components of its 
covariance matrix u vanish except for the momentum- 
momentum component, which coincides with the noise 
kernel. 

The form of the solution is rather simple: all initial 
cumulants of the Wigner function undergo damped os- 
cillations (for the underdamped case) while the ther- 
mal covariance starts from a vanishing value and evolves 
to the asymptotic values corresponding to the thermal 
equilibrium state for the system coupled to the environ- 
ment. We will discuss these solutions more thoroughly 
in Sec. IV. 



D. General Correlations 

Using the initial-value solution of the Langevin equa- 
tion given by Eq. (23) and following the same approach 
as in Ref. [ ], it is straightforward to calculate quan- 
tum correlations between system observables at different 
times. For instance, the symmetrized two-point quantum 
correlation function for position and momentum opera- 
tors in the Heisenberg representation is given by: 



(z(t 1 )z T (t 2 )+z(t 2 )z T (< 1 )) = 
(z(t 1 )z T (f 2 )+z(< 2 )z T (t 1 )) 



(38) 



which with our solutions in Eq. (23) and some basic 
properties of the stochastic Gaussian source, namely 

= and (£(*)£( T )}£ = f(*> r )> wil1 Produce the 
two-time correlation 

((z(< 1 )z T (< 2 )) € ) z =*(< 1 ) < T * T (< 2 ) + ( T T (t 1 ,t 2 ), (39) 



in terms of the two-time thermal covariance 



B. Derivation of the Master Equation 



Jo Jo 

The result for the coincidence-time limit, t\ = t% = 
t, agrees with that of our master equation solution, 
Eqs. (35)-(36), as discussed in Sec. IV A 1. Higher-order 
correlations can be calculated in a similar manner, but 
we can see from the form of our solution in Eq. (35) 
and the Gaussian character of the stochastic source and 
its vanishing mean that only the homogeneous part of 
the solution contributes to cumulants different from the 
second-order one, which are therefore entirely character- 
ized by the initial state of system and the homogeneous 
solutions of the Langevin equation. 



III. MASTER EQUATION 

A. General Theory 

Given the microscopic QBM model of Sec. II A, the 
HPZ master equation for the reduced density matrix op- 
erator p r and for the reduced Wigner function are given 
respectively by 



d_ 
Of 



Pr = - l [H R , p T ] -iT[x, {p, p T }] 

- MD PP [x, [x, p r ]] - D xp [x, [p, p T ]] 



-W r = {H R ,W r } + 2r — (pW r ) 



(41) 



(42) 



+ MD. 



d 2 



tf- 



pp dp 2 1 xp dxdp 



where H R corresponds to the system Hamiltonian with 
f2 2 replaced by a time-dependent frequency fl R (t) ~ £1 2 
whose detailed form, together with that of the time- 
dependent dissipation coefficient T(t) and the diffusion 
coefficients D xp (t) and D pp (t), can be found in Ref. [17]. 

However, as discussed in D, previous derivations of this 
master equation missed a mathematical subtlety concern- 
ing the Green functions of integro-differential equations, 
which renders the existing results for the master equa- 
tion coefficients invalid whenever the nonlocal aspects of 
dissipation become important. In the next subsection we 
provide a compact rederivation of the master equation 
where this issue is properly dealt with, and obtain the 
correct expressions for the coefficients in the general case 
(including the case of nonlocal dissipation). In addition, 
in Sec. Ill C we will provide an analytic expression for 
the solutions of the master equation and show its equiv- 
alence with the result for the state evolution obtained in 
the previous section using the Langevin equation. 



At this point, the quickest derivation of the QBM mas- 
ter equation would merely consist of taking the time 
derivative of Eq. (35) and calculating the inverse Fourier 
transform. Nevertheless, in order to point out the differ- 
ences with previous derivations, which missed the sub- 
tleties of propagators associated with integro-differential 
equations, we will now provide a more traditional deriva- 
tion involving the propagator associated with final- value 
boundary conditions and show that, when done correctly, 
the two are equivalent. We will follow the derivation by 
Calzetta, Roura and Verdaguer (CRV) [19, 20] adapting 
it to our compact notation in terms of phase-space vec- 
tors and matrices. 

We start by considering the stochastic representation 
of the Wigner function 



W T {z,t) = {{8(z(t)-z)) (L 
and differentiate with respect to time: 

d 



Of 



W T (z,t) 



s(t)<*(z(i)-z)>. 



(43) 



(44) 



One can then use the Langevin equation z + H*z = £to 
substitute z(t) and rewrite Eq. (44) as 



d_ 
at 



W r (z,t) 



(45) 



'| f drH(f,T)z(r)^(i)) <5(z(i)-z)^ 



Next, using Eq. (27) one can express z(r) in terms of the 
final value z(t) — z and the propagator <&f(r, r') given 
by Eq. (28). As already pointed out in Sec. II B and 
discussed in detail in D, <frf(r, r') will only be a truly 
advanced propagator [with 3>f (t, t') = for r > r'] when 
considering a strictly local damping kernel, contrary to 
what had been previously assumed. After using Eq. (27) 
we are left with a homogeneous term and two more terms 
involving the stochastic source: 



d_ 
dt 



W(z,t) = V Z T / drH(t,r)*(r,t)zW(z,i) 
Jo 

V Z T Uj^dT J*dr' 'H(t,r) * f (r,r') £(r') S(z(t) -z) 



(46) 



The expectation value of the terms proportional to the 
stochastic source £ can be evaluated with the help of 
Novikov's formula 



<£(r')<5(z(t)-z)) e = 



dr"u{r',T") 



(47) 



Sz(t) 



-iT 



V z 5(z(t)-z) 
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which can be derived by using Eq. (31) and function- 
ally integrating by parts with respect to £. The func- 
tional Jacobian matrix appearing in Eq. (47) can be eas- 
ily obtained by functionally differentiating with respect 
to £(t") the solution of the Langevin equation as given 
by Eq. (23), and one gets 



Sz(t) 



*(t-r"). 



(48) 



Putting these elements together we finally get the follow- 
ing result for the master equation: 

~W T (z, t) = {V Z T n(t) z + V Z T D(t) V z } W T (z, t), (49) 



with the time-local pseudo-Hamiltonian and diffusion 
matrices given respectively by 

U(t) = f drH(i,r)#(T,t), (50) 
Jo 

D(t) = Sy f dru(t,T) * T (t-r)- (51) 
Jo 

Sy j dr f dr'f dr" H(t, r) * f (r, t') i/(t',t") * T (t-r") , 
■/o Jo JO 

and where <frf (r, r') was defined in Eq. (28), and only the 
symmetric part, Sy(M) = (M + M T ) /2, of the diffusion 
matrix contributes to the master equation. These matri- 
ces relate to the conventional representation as follows: 



H(t) = 
D(i) = 



M 



Mfl^(t) 2T(t) 



' 2 Qcp. (^) 



\D xp (t) MD pp (t) 



(52) 
(53) 



The result for the master equation coefficients is ex- 
pressed here in a form analogous to that of previous 
derivations, but this is not the simplest representation. 
We will next proceed to simplify them by eliminating 
the explicit dependence on the time-nonlocal pseudo- 
Hamiltonian H(i, r). 



1. Simplification of the Master Equation Coefficients 

Let us start with the pseudo-Hamiltonian matrix 

H(t) = (H-m)*- 1 (t). (54) 

Taking into account that 3? satisfies the integro- 
differential equation &(t) = — (H • 4 > )(t), the pseudo- 
Hamiltonian can be rewritten as 



nit) = -${t)$-\t). 



(55) 



This new expression for H(t) immediately reveals that 
the homogenous solutions of the nonlocal Langevin equa- 
tion can be equivalently related to the solutions of linear 



ordinary differential equation with time-dependent coef- 
ficients. Indeed, the nonlocal propagator also satisfies 
the dual local equation 



#(t)+K(t)*(t) = 0. 



(56) 



Hence, for local dissipation one would simply have a time- 
independent H and 4>(t) = e^ t ' H , whereas for nonlo- 
cal dissipation H{t) would be time-dependent and <&(t) 
would be given by a time-ordered exponential. 

One can proceed analogously for the diffusion matrix. 
In order to do so we need to simplify the following inte- 
gral: 

J drH(t,T)^(r-T')e(r-r') = J drH(t-r) *(t-t'), 

(57) 



which reduces to 



/t rt — T 

drH(t-T)$(r-r')= / dr H(t-r'-r) #(r) 

= -*(t-r'). (58) 

where we made use of the stationary property of the dis- 
sipation kernel and introduced a simple change of vari- 
ables. Using Eqs. (55) and (58), Eq. (51) can be simpli- 
fied to the following form, which involves terms with at 
most two time integrals: 



D(t) = Sy [ drv{t,T) * T (i-r)- 
Jo 



(59) 



Sy f dr f dr' 
Jo Jo 



d 
dt 



H(t) 



*(i-r) \u{t,t') * t (£-t'), 



where one can clearly see that the second term on the 
right-hand side vanishes for local dissipation, when the 
transition matrix is the exponential matrix e _ . How- 
ever, it can play a crucial role whenever the effects of 
nonlocal dissipation are important, as in the example of 
a sub-ohmic environment of Sec. VI A. 

From our new expression (59) one can see that the 
diffusion matrix can be easily related to the thermal co- 
variance, as given by Eq. (36), and its time derivative. 
Our simplified representation of the master equation is 
then 

~W T (z, t) = { V Z T U(t) z + V Z T D(t) V z } W T (z, t), (60) 

n{t) = -*(t)*- 1 {t), (6i) 

D(t) = \ [U{t) rr T (t) + <r T (t) H T (t) + & T (t)} , 

(62) 

with the phase-space propagator <&(i) given by Eq. (24) 
and the thermal covariance <tt(£) given by Eq. (36). This 
representation contains fewer integrals than the conven- 
tional representation and is completely determined in 
terms of <fr(t) and the noise kernel. 



10 



C. Master Equation Solutions 

In this section we will show that the master equation it- 
self can be solved to produce the same solution as derived 
in Sec. II C. We consider the general master equation 



d_ 

at 



W r = (V Z T D(t) V z + V Z T U{t) z) W t . 



(63) 



This is a hyperbolic second-order partial differential 
equation (PDE). The equation is not separable in time 
nor phase-space. The nature of the PDE suggests taking 
a Fourier transform of the phase-space variables as the 
derivatives are of higher order than the algebraic param- 
eters. Furthermore, not only does a Fourier transform 
reduce the PDE to first order, but the computation of 
expectation values also becomes trivial since we are then 
working with the characteristic function of the distribu- 
tion. 

The Fourier transform is defined as 



/+oo p+oo 
dx dp 
■oo J—oo 



and it exhibits the usual properties: 



7(z), (64) 



+00 /> + OG 

dx dpq r *f{z) 

— oo J—oo 



The master equation becomes then 

I- + ^-HM, ) Wr = k T D kW r . 
at I 



(65) 



(66) 



where W r = ^{M^} and the normalization of W T (z, t) 
implies W T (0,t) = 1. 

From Eq. (66) it is clear that if the master equation 
coefficients asymptote to constant values, then we will 
have a stationary Gaussian solution in the late-time limit 
given by 



(67) 



with cr^f uniquely determined by the Lyapunov equation 



U —OO I _.oo /T/T 

Ttoo cr T ~i~ °T n oo 



2D r 



(68) 



To zeroth-order in the system-environment coupling, this 
corresponds to the free thermal state of the system. It 
is also reasonable to believe that more generally this cor- 
responds to the thermal state of our system coupled to 
the environment (i.e. the reduced density matrix of the 
thermal state of the whole system including the system- 
environment interaction). For arbitrary systems this has 
been proven to second order in the system-environment 
coupling (here first order in damping, e.g. 70) [? ]. 



1. Method of Characteristic Curves 

The method of characteristic curves involves looking 
for parameterized curves in the domain (t, k) along which 



the first order PDE becomes a set of first-order ODEs. 
For each one of those curves we have 



W r [k,t]=W r [k(T),t(T)], 



W x 



^Wr 
dT dt * 



dk l 



V k W r , 



(69) 
(70) 



dr dr dt dr 

Next, we attempt to match the right-hand side of Eq. (70) 
to the left-hand side of Eq. (66). This results in a system 
of ODEs in the parameter r. We will look for curves 
that synchronize with the initial time so that t(Q) = 0, 
k(0) = ko. The solution for the parameterization of the 
time coordinate is simple: 



dt 
rf7 



= 1 => t(r) = r. 



(71) 



On the other hand, finding the parameterization for the 
Fourier transform of the phase-space variables is a bit 
more involved. It is characterized by the linear ODE 
system 



dr 



k T (r) = +k T (r)^W 



and its solutions can be written as 
k(r) = * fe (r)k , 



(72) 



(73) 



where 3>fc (r) is the matrix propagator associated with the 
transpose of Eq. (72) and equals the identity matrix at 
t = 0. For local dissipation, H is time independent and 
the propagator is simply given by *&J(t) = e +rH , which 
equals $ _1 (r). Such a relation between the matrix prop- 
agator of the integro-differential Langevin equation (19) 
and the local equation (72) actually holds in general. In- 
deed, taking into account Eq. (55), it follows that the 
propagator for the characteristic curves (r) must sat- 
isfy the equation 

| : *J(r) = -*J(r)*(t)*- 1 (t), (74) 

which is equivalent to the relation 
d 

Jr 

Together with 4>£(0) 3>(0) = I, since both *fe(r) and 
3>(t) equal the identity matrix at the initial time, this 
implies that *J(r) = * _1 (r). 

We now have the rules for transforming back and forth 
between the domain coordinates (t, k) and the character- 
istic curve coordinates (r, ko); ko uniquely specifies each 
characteristic curve parameterized by r. Using these re- 
sults, we can immediately apply the method of charac- 
teristic curves to solving Eq. (66) as follows: 



(#£(t)S(t)) =0. 



(75) 



dr 



WMr),t(r)} = -k T D(t)kW r [k(r),((r)], 



(76) 



dr 



W r [* fc (r)ko,t(T)] = -k£#j£(T)D(T)* fc (r)ko 

x W r [* fc (r)k ,t(r)]. (77) 
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The last equation is a linear ODE whose solution can be 
easily found to be 

W r [* fe (r)k 0! r] = (78) 

W r [k ,0]e-- f » Tdr '( k «*' (T ' )D(r ' ) * t(T ' )k|, ) I 

where W r [ko,0] is the initial characteristic function at 
t = 0. We can now express the solution back in terms of 
k and $ to get the final result 



Wr[M] = W r 



* T (i) k, 



,-ikVr(t)k 



with thermal covariancc defined 



(T T (t) = 2 [ dT*(i,T)D(r)* T (t,T), 
Jo 



(79) 



(80) 



and note that 3>(t, r) here does not have time- 
translational invariance for nonlocal dissipation, where 

sec the discussion in 

D. 



2. Equivalence with the Result from the Langevin Equation 

We have shown that the form of the solution from 
the master equation is equivalent to that derived from 
the Langevin equation in Sec. II B. What remains to be 
shown is that the thermal covariances are indeed equiv- 
alent. To do this one can differentiate Eq. (80) with 
respect to time and get the following result: 

& T (t) = -H{t) a T (t) - tr T {t) H T {t) + 2 D(t). (81) 

This equation is also satisfied by the thermal covariancc 
expression directly derived from the Langevin equation, 
as can be seen from Eq. (62). Furthermore, the ther- 
mal covariances given by Eqs. (80) and (36) both have 
vanishing initial conditions: <x T (0) = 0. Therefore, since 
they are both solutions of the same ordinary differential 
equation and have the same initial conditions, they must 
be equivalent. 



where <&(i) is the phase-space propagator for the 
Langevin equation defined in Eq. (22). 

The solution in Eq. (82) consists of two factors. The 
first one tends to unity in the long time limit and en- 
codes the disappearance of the initial state (we will call 
it the death factor). The second factor describes the ap- 
pearance of a Gaussian state that evolves in time and 
tends asymptotically to a state that corresponds to ther- 
mal equilibrium (we will refer to this as the birth factor) . 
Assuming dissipation, all initial distributions evolve to- 
wards this final Gaussian state, with thermal covariance 
(T-rif)- This state does not look like the thermal state 
of a free harmonic oscillator because of the coupling to 
the environment. It more likely results from considering 
the thermal equilibrium state for the whole system (sys- 
tem plus environment) including the system-environment 
interaction, which gives rise to a non-trivial correlation 
between them, and tracing out the environment. 

The death factor contains the information on the ini- 
tial conditions; it describes the gradual disappearance 
of the initial distribution and it is always temperature 
independent. The free evolution of the Wigner func- 
tion corresponds to rotation in phase space (when prop- 
erly rescaled) at constant angular velocity. Dissipation 
will modify this rotation to inspiralling of the trajecto- 
ries down to the origin, or decay to the origin without 
completing a full rotation in the case of overdamping. 
More generally, for nonlocal dissipation the trajectories 
will correspond to those of a parametrically damped os- 
cillator, which in some cases could be quite complicated. 

The birth factor describes the complicated birth and 
settlement of a state of thermal equilibrium. This fac- 
tor is always Gaussian with a covariance matrix given by 
Eq. (83) , which involves a convolution of the noise kernel 
with propagators that reflect the natural oscillatory de- 
cay of the system. This covariance matrix vanishes at the 
initial time and tends at late times to an equilibrium co- 
variance matrix which can be easily determined from the 
Lyapunov equation (68). The thermal covariance matrix 
is always positive definite. 



IV. EVOLUTION OF STATES 



Trajectories of the Cumulants 



A. General Solutions 

Whether derived via the Langevin equation in Sec. II C 
or solving the master equation in Sec. Ill C, the evolution 
of the system state is most easily represented in terms of 
the characteristic function (the Fourier transform) of the 
reduced Wigner distribution: 



W r [M] = W r [* T (t)k,oj e -5k T ^T(t)k 
with the thermal covariance <xr(i) given by 



(82) 



(T T (t)= ( dT [ dT' ®(t-T)v(T,T')$ T {t~T'), (83) 

Jo Jo 



As we have already mentioned, the Fourier transfom 
of the reduced Wigner function corresponds to its char- 
acteristic function, from which the correlation functions 
for the phase-space variables can be easily derived using 
Eq. (65). The general expressions for the cumulants can 
be obtained straightforwardly from the logarithm of the 
reduced Wigner function in Fourier space as follows: 



oo n 



(84) 



n=l 



where k H denotes the components of the vector k and we 
used the Einstein summation convention for pairs of re- 
peated indices (i.e., it is implicitly understood that a sum 
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Si =1 should be preformed over each pair of repeated in- 
dices «W is the n th cumulant and acts as a tensor of 
order n contracted with n copies of k. Using the result 
for W T (t, k) from Eq. (82) we have 



E 



1 n co n 

i=i 

,kV(«)k, (85) 



n=l 

2 



where k^™' 1 - n (0) are the cumulants associated with the 
initial distribution. Eq. (85) implies 

n - ■ 

(*) = (°) n * (* T w) Ji ' + ^2 4 iia (*). 



8=1 



. (86) 

We can see that the only cumulant with a non- vanishing 
asymptotic value, which is a consequence of the thermal 
fluctuations, is the covariance matrix (with n = 2). The 
closely related second momenta of the distribution are 
given by 



(zz T )(t) = *(t) (zz T ) qo * T (t) + <x T (i), 



(87) 



where 



denotes the expectation value with respect 



to the reduced Wigner function at the initial time, 1 as 
defined in Eq. (30). All other cumulants experience what- 
ever oscillatory decay is inherent in the homogeneous so- 
lution of the Langevin equation. In particular, the ex- 
pectation value 



<z)(t) = *(t) (z) 



follows a trajectory like that plotted in Fig. 1 for local 
dissipation, where one can see that the trajectory of the 
expectation values (x) , (p) for any initial distribution in- 
spiral into the origin. This captures the behavior of Gaus- 
sians plotted by Unruh and Zurek [16]. 



2. Thermal Covariance 

As we have seen, the only additional quantity that 
needs to be calculated besides the propagator is the ther- 
mal covariance. Here we discuss the full-time evolution 
of the thermal covariance, which can be most easily ob- 
tained from Eq. (83). Using the addition formula for the 



(p) 

MSI 




FIG. 1. The trajectory of the expectation values (x), (p). 



argument of the cosine function appearing in the defi- 
nition of the noise kernel, one obtains the following ex- 
pressions for the components of the thermal covariance, 
which only involve calculating a single time integral be- 
sides the integral over frequencies: 



of (t) = 



du coth(|^) [G(t) * cos {ujt)f (89) 

+ [ duI(uj)coi\\(—\ [G(t) *sm(ujt)} 2 , 
Jo \2T ) 

(90) 



<4"(t) = -MoT®, 



af{t) 



cko I(ui) coth 



2TJ 



MG(t) * cos (tot) 



dw/(w)coth(|^ MG(t) * sin (tut) 



(91) 



Note that the expectation value of any phase-space function with 
respect to the reduced Wigner function is equivalent to a quan- 
tum expectation value with respect to the corresponding reduced 
density matrix where the arguments x and p of the phase-space 
function are promoted to operators and the Weyl ordering pre- 
scription is employed. In particular, for the second-order cu- 
mulants this corresponds to considering symmetrized two-point 
quantum correlation functions. 



These results are expressed in terms of Laplace convolu- 
tions of the propagator with sinusoidal functions which 
become trivial in Laplace domain, although one must 
eventually transform back to compute the squares. More- 
over, integrating by parts in the Laplace convolutions and 
taking into account that G(0) = 0, the momentum co- 
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variance can be expressed in the alternative form 
a P p (t) = [ (1ujuj 2 I(uj) coth( — ) [MG(t) * cos (ujt)} 2 

.In V 27/ 



dww 2 /(w)coth(— ) [M G(t) *sin(wt)] 

n V27V 
2 ,.((\\ /^</j-\2 



M z v{Q)G{ty 



(92) 



which is completely analogous to that for the position 
covariance, but with an effectively higher-order spectral 
density due to the additional factor of w 2 , plus a simple 
cut-off sensitive transient term which decays with the 
characteristic relaxation rate. It becomes then obvious 
that the momentum covariance will contain the domi- 
nant contribution to any potential ultraviolet sensitivity 
of the thermal covariance, whereas the position covari- 
ance will contain the dominant contribution to any pos- 
sible infrared sensitivity. 

In order to compute the evolution of the thermal co- 
variance, especially when calculating it numerically, it 
is often convenient to use the following alternative ex- 
pressions, which can be derived by differentiating with 
respect to time the xx and pp components of Eq. (83): 

a%?(t) = 2G[t)[v(t)*G(t)}, (93) 
d 



'(t) = 2M 2 G(t)j t {v(t)*G(t)} 
M 



(94) 



= (*) = MG{t) [v{t) * G(t)} , (95) 

where the convolution of the propagator with the noise 
kernel should be performed before the frequency integral 
of the noise kernel. This will typically result in expres- 
sions more amenable to numerics since one can avoid in- 
creasingly oscillatory integrands. 

For odd meromorphic spectral functions the frequency 
integral can be evaluated by contour integration (and the 
residue theorem) using the rational expansion of the hy- 
perbolic cotangent 



coth 



,27V = 



2T 



E 



2irT 



7T <■ — ' h2 

fc=l K 



(96) 



One should then be left with a sum of terms ratio- 
nal in the Laplace domain, which can be contracted 
into digamma or harmonic-number functions [respec- 
tively tp(z) or H(z)], which are asymptotically logarith- 
mic. When transforming back to the time domain, the 
residues of the hyperbolic cotangent additionally give rise 
to products of rational functions of k with e - 27rT * fc . These 
terms contain all effects which decay at temperature- 
dependent rates and can be expressed in terms of Lerch 
transcendent functions, <P(z, 1, e~ 27rTt ) , which are useful 
for numerical calculations but not particularly insightful. 

Fortunately, one can also derive a simple analytic ex- 
pression for the late-time thermal covariance, as shown 
in E: 



erT'(oo) 



dco I(uj) coth 



2T 



|GM| 5 



1 

M 2 w 2 



(97) 



which reduces the calculation of late-time uncertainties 
to a single integral. This relation confirms that for late 
times the momentum covariance has precisely lu 2 more 
frequency sensitivity in its integrand. 



3. Linear Entropy 

In this subsection we investigate the linear entropy [31], 
which can be easily obtained from the Wigner distribu- 
tion as follows: 



S L = 1 - Tr(p 2 ) = 1 - 2tt Jd 2 zW 2 {z,t). 
In Fourier space it becomes 

s L =i-^ yd 2 kiw r (k,t)i 2 , 

and using the result in Eq. (82) we finally get 



St 



1 



2n 



cfk 



(o,* T (t): 



= -k T <T T (t)k 



(98) 



(99) 



(100) 



At the initial time the linear entropy is that of the initial 
state, and at late times it tends to 



St 



1 



1 



2ydet7. 



(101) 



Alternatively, one can express the linear entropy in 
terms of an integral of the Fourier-transformed reduced 
Wigner function at the initial time by introducing the 
change of variables ko = $ T (t) k. Eq. (100) can then be 
written as 



Sr. 



1 



d 2 k 



1 



2tt J ~ " u det[*(f) j 
x e -kT#- 1 (t) o - T (t)#- T (t)k 



|W r (0,k ) 



= 1 - 



1 



2 A /det[ < r T (t)] 



d 2 k |W r (0,k ) 



x n(o, ^ T (t) cr T \t) #(t);ko) , (102) 

where N(fi, cr; ko) is a normalized Gaussian distribution 
for the variable ko with mean \i and covariance er. For 
small times this integral is similar to that for the initial 
state, whereas for long times the normalized Gaussian 
distribution becomes increasingly close to a delta func- 
tion. 

For a Gaussian initial state 

W r (0, k ) = exp (-kfo-o k - i k r (z) ) , (103) 
the integral in Eq. (100) can be explicitly computed: 



Sr 



2tt 



d 2 ke" kT (* (t)CTo * T(t)+CTT(t) ) k 



= 1 - 



2 A /det *(t)er * T (i) + cr T (t) 



(104) 
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For these Gaussian states, reasonable linear entropy is 
synonymous with reasonable uncertainty functions (i.e., 
the linear entropy will be positive if and only if the 
Heisenberg uncertainty principle is satisfied). We will 
find that the late time uncertainty is well behaved. The 
uncertainty at the initial and intermediate times should 
not violate the Heisenberg uncertainty principle either. 



4- Decoherence of a Quantum Superposition 

In this section we will illustrate how one can get 
a useful qualitative picture of the phenomenon of 
environment-induced decoherence from the the solutions 
of the master equation given by Eqs. (82)-(83). In or- 
der to do that we will consider a quantum superposition, 
\tp) = + \ip~))/VK, of a pair of states \ip±) which 

correspond to a pair of Gaussian wavefunctions in posi- 
tion space separated by a distance 28 x and where K is 
an appropriate normalization constant. Specifically, we 
have 



ij>±(x) = ip {xT$x), 



ik(x) = yjN(0,(rS*;x). 



(105) 
(106) 



where N(ji, er 2 ; x) is a normalized Gaussian distribution 
for the variable x with mean /i and variance c 2 , and i/jq(x) 
is a reference Gaussian state centered at the origin. 

Taking into account the definition of the Wigner func- 
tion, 



W(x,p) = 



1 

2^ 



dye™p(x-y/2,x+y/2), (107) 



and applying it to the density matrix p(x,x') 
(x\ip) (ip\x r ) we get 

W(z) = — \w+(z) + WL(z) + 2cos(24p) W (z) 



K 



(108) 



where W + , WL and W are respectively the Wigner func- 
tions of the states |V>-) and |-0o). This Wigner 
function, plotted in Fig. 2, exhibits oscillations of size 
1 /6 X along the p direction. These oscillations are closely 
connected to the coherence of the quantum superposition 
(and the existence of non-diagonal terms in the density 
matrix) and are absent in the Wigner function for the 
incoherent mixture W(z) = (l/2)[W + (z) + W-(z)]. 

In this context the decoherence effect due to the inter- 
action with the environment corresponds to the washing- 
out of the oscillations in the reduced Wigner function as 
it evolves according to the master equation. This can 
be seen rather simply from the result for the solutions of 
the master equation obtained in this section and given 
by Eqs. (82)-(83). Taking into account that the inverse 
Fourier transform of Eq. (82) corresponds to a convolu- 
tion of the homogeneously evolving initial state and a 
Gaussian function with the thermal covariance (Txit) as 




FIG. 2. Wigner function associated with a state \ip) = 
(IV'i) + \ip2})/V~K which corresponds to the coherent quan- 
tum superposition of two Gaussian wavefunctions in position 
space shifted by a distance 8 X . 



its covariance matrix, the Wigner function can then be 
expressed as 



W T (t,z) = / dz 



sN{ 0,(T T (t);z-z') 
det[*(t)] 



W r (Q^-\t)z') 



(109) 

where the thermal Gaussian acts as a Gaussian smear- 
ing function which starts as a delta function at the ini- 
tial time and broadens with the passage of time until 
it eventually reaches its asymptotic thermal-equilibrium 
value. Therefore, several aspects will be at play. On 
the one hand, the initial state evolves as a phase-space 
distribution with trajectories corresponding to the ho- 
mogeneous solutions of the Langevin equation (19) and 
with the same qualitative behavior depicted in Fig. 1 
for the trajectories of (x) and (p). On the other hand, 
by diagonalizing <Tt(<) at each instant of time one gets 
the principal directions and the widths ((71,02) of the 
Gaussian smearing function, which will average out any 
details of those sizes along the corresponding directions. 
When <Tt(*) along the direction of the interference os- 
cillations of the Wigner function becomes comparable to 
their wavelength, they get washed out and the Wigner 
function becomes equivalent to that of the completely 
incoherent mixture. The time it takes for this to happen 
is known as the decoherence time tdec- 

Knowledge of the qualitative behavior of <Tr(f), com- 
bined with the fact that the phase-space distribution 
det[€>(i)] _1 W T {0, ^ rl (t) z') is rotating with the charac- 
teristic oscillation frequency and shrinking with the char- 
acteristic relaxation time is all that one needs to un- 
derstand how different initial states decohere as time 
goes by. In particular, if the decoherence time-scale, 
given by tdec, is much shorter than the characteristic 
oscillation period and the relaxation time (but suffi- 



15 



ciently longer than 1/A), one can approximate the phase- 
space distribution by the initial reduced Wigner func- 
tion (after any possible initial kick). For instance, for 
an Ohmic environment in the high-temperature regime 
one can, under those circumstances, approximately take 



.pp/ 



'rp it) ~ D™t with ~ 2Mj T and from the con- 
dition a/ Opp (t) ~ 1/S X obtain an estimated decoherence 
time idee ~ ^-/(^MjoTS^.), in agreement with the stan- 
dard result for this situation [32, 33]. On the other hand, 
if M, 70 or 5 X are very small t& cc can become compara- 
ble or larger than the dynamical timescales 1/17 or 1/7, 
and the previous estimate can no longer be applied be- 
cause one needs to take into account the evolution of 
<Ty(t), which is then less simple (it will roughly oscillate 
with frequency 17 around a central value which increases 
with a characteristic timescale 1/7 until it approaches 
the asymptotic thermal value), as well as the rotation 
and shrinking of the initial Wigner function under the 
homogeneous evolution. Note also that if we had consid- 
ered an initial superposition of Gaussian states peaked at 
the same location but with different momenta, which cor- 
responds to a Wigner function along the position rather 
than momentum direction, the decoherence time would 
typically be much longer, since a^ x (t) vanishes at the 
initial time and grows with a characteristic timescale of 
order 1/17. In that case, the rotation of the Wigner func- 
tion becomes important since the oscillations can then 
be averaged out due to the larger values of a^(t). 

The zero-temperature regime for an Ohmic environ- 
ment is also qualitatively different. There is a substan- 
tial contribution to 0jf (t) from a jolt of the diffusion 
coefficient D pp for times of order 1/A. However, this is 
actually regarded as an unphysical consequence of hav- 
ing considered a completely uncorrelated initial state for 
the system plus environment, and this kind of highly 
cut-off sensitive features at early times of order 1/A 
should disappear if one considers a finite (cut-off inde- 
pendent) preparation time for the initial state of the sys- 
tem coupled to the environment [34]. For further discus- 
sion on this point as well as a possible way of avoiding 
these spurious effects and generating a properly corre- 
lated initial state by using a finite switch-on time for the 
system-environment interaction see C 2. For sufficiently 
weak coupling, M, or 5 X , t^ cc can become comparable 
or larger than the relaxation time more easily than at 
high temperatures since the components arxit) are much 
smaller in this case. For example, the asymptotic ther- 
mal value of Opp is of order MCI (for weak coupling), 
much smaller than the high-temperature results, which 
is of order AIT. In such situations, the main effect of con- 
sidering a sufficiently long time is through the shrinking 
of det[*(t)] _1 W T (0, * _1 (i) z') and the size of its oscilla- 
tions. 

We have focused in this subsection on describing the 
qualitative features of the environment-induced decoher- 
ence of an initial coherent superposition that can be eas- 
ily inferred from our general result for the evolution of 
the reduced Wigner function. A much more quantitative 



study is possible by using the exact analytical results 
for the diffusion coefficients and, especially, <xt(£), which 
will be presented in Sees. V and VI. We expect agreement 
with the numerical results obtained in Ref. [33], although 
significant deviations may appear when the nonlocal ef- 
fects of dissipation are important (such as in the sub- 
ohmic case) since previously obtained master equations 
are not valid in those regimes. 



B. Late-Time Dynamics 

We now focus our attention on the dynamics generated 
by the stationary limit of the master equation, assuming 
that one exists. For an Ohmic spectrum with a large cut- 
off the pseudo-Hamiltonian H will reach its asymptotic 
value within the cut-off timescale, whereas the diffusion 
D within the typical system timescales (although certain 
terms contributing to the diffusion coefficients will decay 
at a temperature-dependent rate whenever this is faster) ; 
see Sec. V for a detailed analysis of all these questions. In 
the weak-coupling regime this leaves the majority of the 
system evolution within this late-time regime wherein the 
master equation is effectively stationary. However, the 
existence of such a regime is not guaranteed in general. 
For instance, in the sub-ohmic case the evolution can be 
persistently nonlocal and the effectively local late-time 
regime discussed here need not exist, as will be shown in 
Sec. VIA. 



1. Late-Time Propagator 

If the late-time stationary limit of the master equa- 
tion exists, the late-time pseudo-Hamiltonian operator 
will take the form 



n 







" M 

M17 R 2r 



(110) 



and can be effectively represented as arising from the 
propagator 



G R (s) 



M 



s 2 + 2Ts + 17 R ' 



-ri 



(111) 
(112) 



with 17r = \/l7 R — r 2 . This effective propagator Gr(£) 
is not equivalent to the late time limit of the true prop- 
agator G(t), but they should share the same asymptotic 
dynamics. Specifically if one can take the asymptotic 
expansion 



G(t) = G 00 (t) + 5G(t), 



(113) 



where Goo(i) contains the asymptotic limiting behavior 
and SG(t) contains the early time corrections, which de- 
cay faster at late times, then Goo {t) should directly yield 
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£7r and r in its arguments, although a phase and am- 
plitude difference between G 00 (t) and Gr(£) may exist. 
This can be rigorously justified if 7 (s) and, thus, G(s) 
are rational, which implies that the time dependence 
of G(t) corresponds to damped oscillations with various 
timescales. On the other hand, the sub-ohmic spectral 
distribution that will be studied in Sec. VI A provides a 
pertinent counter-example [in that case G(t) decays as 
a negative power-law rather than exponentially] which 
shows that this situations does not necessarily exist when 
the spectral density function is not meromorphic. 

If we indeed have a rational spectral density, then from 
the nonlocal propagator one only needs to solve the char- 
acteristic equation 

/ 2 + 2 7 (/) / + f! 2 = 0, (114) 

to obtain all the rates / associated with the propagator 
(this is the same equation whose roots need to be found 
when decomposing the propagator in Laplace domain 
into simple fractions). From Eq. (114) and the positivity 
of the damping kernel, it follows that the real part of / 
will always be negative definite. Those with the smallest 
real part in absolute value give the late-time coefficients: 
the real part corresponds to — T and the imaginary one to 
Qr. A specific example can be found in Sec. V, where the 
Ohmic case with a finite cut-off is studied in detail. On 
the other hand, if one treats the system-environment in- 
teraction perturbatively, one can show that the late-time 
weak-coupling coefficients take the following form: 



/± = -r±tn R , 

r = Re[7^)]+0( 7 2 ), 
il R = n-Im[7(ifi)] +C(7 2 ), 



(115) 
(116) 
(117) 

which is in agreement with the results for the weak- 
coupling master equation obtained in Ref. [? ]. Any 
additional timescales would then be perturbations of the 
cut-off or other timescales intrinsic to the spectral func- 
tion. 

It should be noted that in general the late-time prop- 
agator discussed here cannot be employed to calculate 
the diffusion coefficients or the thermal covariance, not 
even at late times. This is because both quantities eval- 
uated at an arbitrary time t get non-negligible contribu- 
tions involving the propagator at early times, as can be 
seen for instance from Eqs. (59) and (83). Nevertheless, 
one can still employ the late-time propagator to obtain 
the late-time evolution of the thermal covariance (and 
the diffusion coefficients) provided that one already has 
an accurate result for its constant asymptotic value [ob- 
tained for example with Eq. (97)], as will be illustrated 
next. In addition, one can also use the propagator Gr(£) 
given by Eq. (112), which corresponds to the limit of lo- 
cal dissipation, to calculate the thermal covariance and 
diffusion coefficients for an Ohmic environment with a 
sufficiently large cut-off, since in that case the contribu- 
tion from the extra early-time term of the propagator can 
be neglected when calculating these quantities for times 
later than A -1 , as will be shown in Sec. V. 



2. Late-Time Diffusion and Covariance 

Given late-time master equation coefficients which 
have all taken their asymptotic values, one can show that 
the evolution of the covariance in that regime is given by 

cr(t) = o-5? + *(t-*0 [<r(ti) - er^?] # T (i-ii), (118) 

which is a solution of Eq. (81) as long as one assumes 
H{t) and D(t) to be time-independent after some time ij 
in the late-time regime. Note that we have assumed that 
the master equation coefficients reached their asymptotic 
values much faster than the relaxation time (as illus- 
trated in F with the example of the ohmic distribution, 
this may be the case for finite temperature, but not nec- 
essarily so for zero temperature). 

The asymptotic value of the late-time thermal covari- 
ance <r^? has been reduced to a single integral in E. From 
this single integral formulation, it is actually easier to 
obtain first c^?, and then obtain the late-time diffusion 
coefficients using the Lyapunov equation (68). However, 
it is interesting to note the inverse relation 
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(119) 



for the following reason. As we have pointed out in 
Sec. IV A 2, only the momentum covariance can contain 
the highest frequency sensitivities. From the Lyapunov 
solution we can sec that the regular diffusion coefficient 
would also contain such high-frequency sensitivities as 
it alone determines the late-time momentum covariance. 
Therefore, the anomalous diffusion coefficients must act 
as an "anti-diffusion" coefficient in keeping the position 
covariance free of such sensitivities. On the other hand, 
only the position covariance can contain the lowest fre- 
quency sensitivities and these must, therefore, be entirely 
contained in the anomalous diffusion coefficient if they 
exist. 

In summary, any specific features of the initial dis- 
tribution decay away and at late times the state tends 
generically to a Gaussian with a covariance matrix given 
by Eq. (97). As follows from Eq. (87), the late-time posi- 
tion and momentum uncertainties are, therefore, entirely 
given by the asymptotic values of the thermal covariance: 

(120) 
(121) 



(Ax) 2 = (<r~) a 
(Ap) 2 = („§?), 



V. OHMIC CASE WITH FINITE CUT-OFF 

A. The Nonlocal Propagator 

The arguably simplest example of ohmic dissipation 
with finite cut-off that one can construct corresponds to 
the following damping kernel: 



7 (s) 



7o 



1 + X 



(122) 
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This damping kernel is constant at frequencies much 
smaller than the cut-off, but vanishes in the high fre- 
quency limit. The corresponding spectral density also ex- 
hibits a rational cut-off function, which decays quadrat- 
ically for large frequencies: 



I {id) = — M70 id 
7T 



fid 

1+ (a 



(123) 



Calculating the Green function amounts to factoring 
a cubic polynomial. Specifically, one needs to factor 
(s 2 + Sl 2 )(,s + A) +270 As in the denominator of the Green 
function G(s). For the underdamped system the effect of 
a large finite cut-off is to shift the system relaxation and 
oscillation timescales slightly: 



7* = 7o 

n 2 = 



A 



A - 27* 



n 2 



(124) 
(125) 



and to add an additional relaxation timescale comparable 
to the cut-off: 



A. = A - 2 7 ,. 



(126) 



If we parametrize everything in terms of these phe- 
nomenological frequencies, the Green function for the 
fully nonlocal damping kernel can always be expressed 
as 



G{s) 



1 



A 



M {s + AO (s 2 + 2 7 „s + n 2 ) ■ 



(127) 



without the need to explicitly factor a cubic polynomial, 
while the original parameters are given by 



7o 

n 2 

A 



a; + 2 7 a + n; 

(A, + 2 7 *) 2 
A. o2 



-7* 



A* + 2 7 * 
A* + 2 7i . 



(128) 

(129) 
(130) 



then we never have to actually factor the cubic polyno- 
mial. 

After using partial fraction decomposition in Eq. (127), 
one can easily transform back to the time domain and 
obtain the exact propagator for the nonlocal case: 



G(t) = 



A 2 + n 2 



(A*- 7 *)2 + fi2 



G R (t)- 



27, 



Al+n 2 



Gr(*)- 



,-A*i 
M 

(131) 



where Gn(t) is the late-time local propagator introduced 
in Eq. (112). Note that as long as A* > 7* the term 
proportional to e~ A ** can be neglected at sufficiently 
late times, when the terms involving Gr(£) dominate. 
This corresponds to the late-time regime discussed in 
Sec. IV B 1 [the term proportional to Gn(t) simply causes 



a phase shift] and the late-time master equation coeffi- 
cients are, therefore, 



r = 7*, n R = n*. 



(132) 



In the high cut-off limit one recovers the usual coefficients 
70 and Q. Furthermore, in that limit one can approxi- 
mate G(t) by Gn(t) since the extra terms are suppressed 
by inverse powers of A 2 . For G(t) this is true even at 
arbitrarily early times of order A -1 : although the ex- 
ponential factor is not suppressed, the prefactor 1/A 2 
is sufficient to suppress its contribution to G(t). This 
is not true, however, for G(t) (or higher-order deriva- 
tives), which also appears in 4>(t). From Eqs. (36) and 
(24) we can see that the component involving G(t) does 
not contribute to the thermal covariance, but whether it 
contributes to its time derivative &r(t) as well as to the 
diffusion coefficients, which are related to <tt(£) through 
Eq. (62), is a bit more subtle. In order to analyze this 
point it is convenient to consider Eq. (59). On the one 
hand, the time derivative acting on <f>(t— r) in the second 
term on the right-hand side of that equation will give 
rise to G{t — r) and an unsuppressed contribution from 
e -A,(t-T)^ [Analogously to what was explained above for 
Eq. (36), there is no contribution from the components 
of the transition matrix involving G(t), and it can only 
arise when time derivatives act on other components.] 
On the other hand, the additional time integral in that 
term when considering such a contribution will generate 
an extra 1/A* factor as compared to the first term on the 
right-hand side of Eq. (59). Thus, the final conclusion is 
that we can use the approximate local propagator G^ (t) 
to calculate the diffusion coefficients at arbitrary times in 
the large cut-off limit. Comparison of the results evalu- 
ated using the exact expressions and plotted in Sec. VB 
and the approximate results for the large cut-off limit 
also support this conclusion. 

We close this subsection with a brief discussion of the 
possible dissipative regimes when considering finite val- 
ues of the cut-off in our spectral function, since the pres- 
ence of this new scale can give rise to a richer set of 
possibilities. For our rational cut-off function we have 
three different dissipative regimes corresponding to the 
three shaded regions in Fig. 3. The boundary between 
different regions corresponds to the values of the param- 
eters for which a pair of roots of the denominator of G(s) 
degenerate and change character, i.e. they change from 
a complex conjugate pair to two real ones. Atop the di- 
agram where A 3> f2, lies the regime of local dissipation, 
whereas along the bottom of the diagram where A <C £7, 
lies an effectively sub-ohmic regime as A becomes an IR 
cut-off. The white shaded vertical stripe to the left lies 
completely in the weak coupling regime and constitutes 
the underdamped regime. This regime is as described 
previously with slowly decaying oscillations and a cut- 
off-dependent decay rate. The grey shaded middle region 
denotes the overdamped regime. This regime is also anal- 
ogous to that of the simple and overdamped harmonic 
oscillator but with an additional cut-off-dependent decay 
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FIG. 3. Dissipative phases for Ohmic damping with finite 
rational cut-off. From left to right they are underdamped in 
white, overdamped in grey, and strong coupling in black. 



rate. The black shaded region to the right denotes a new 
nonlocal strong- coupling regime that emerges for a suffi- 
ciently strong coupling (such that 70 is large compared 
to the cut-off). Specifically, as derived from Eqs. (128)- 
(130), the relevant scales for this regime in the limit of 
very strong coupling are 



A* = 



7* = 



27o 
A 

2 ~ 



4A 7o 2 



O 



47o 



o 



n+ = 2A 7o + n 2 + o 



(133) 
(134) 
(135) 



Hence, we can see that one has moderately damped, rapid 
oscillations plus an additional slow decay rate. 



B. Initial Jolts 

Early studies by Unruh and Zurek [16] as well as HPZ 
[17] already revealed that at low temperatures the nor- 
mal diffusion coefficient D pp (t) of an ohmic environment 
exhibited a strong cut-off sensitivity for very early times 
of order 1/A. As shown in the next section and F, in 
the large cut-off limit where the use of the local propa- 
gator is a good approximation one can obtain relatively 
simple analytic results. They confirm that for zero tem- 
perature the normal diffusion coefficient, which vanishes 
at the initial time, exhibits an initial jolt with an ampli- 
tude of order A peaked around a time of order 1/A and 



then decays roughly like 1 ft (for times much earlier than 
1/Q and I/70). 

Alternatively, one can obtain the exact analytic results 
for finite cut-off by computing &r(t) using Eqs. (93)-(95), 
as explained in Sec. IV A 2. The resulting expressions are 
rather lengthy and not particularly insightful, and will 
not be reported here, but they have been employed to 
plot some examples of exact results for the diagonal com- 
ponents of &r(t) and <xt(£) in Figs. 4 and 5. From the 
different components of the thermal covariance and its 
time derivative one can obtain the diffusion coefficients 
using Eq. (81), and in particular one can see from Fig. 5 
the presence of the jolt mentioned above. 





FIG. 4. Exact thermal covariance dynamics for • normalized 
position uncertainty MQ+ax X (t) and ■ • ■ normalized momen- 
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FIG. 5. Same plot as in the previous figure, but with a much 
larger time resolution, which reveals the presence of the initial 
jolt in o- pp {t) peaked around t ~ 1/A*, while a xx (t) and <7 xx (t) 
remain essentially zero at those timescales. 

It is important to emphasize that such kind of behav- 
ior, as well as an associated rapid growth of o~ pp (t) and 
a slower growth of o~ xx (t) (which eventually decays ex- 
ponentially within the relaxation time-scale 1/r) until 
they both reach values which depend logarithmically on 
A for large values of A, is a consequence of having started 
with a completely uncorrelated initial stated. A possible 
way of generating a properly correlated initial state is 
by smoothly switching on the system-environment inter- 
action within a time-scale much longer than 1/A, but 
longer than the other relevant time-scales of the system. 
This is discussed in some detail in C 2. It also contains a 
number of technical details concerning the effects of the 
switch-on function appearing in the dissipation kernel, 
which can be be mainly reabsorbed in redefinition of the 
initial sate. The key point, however, is the role played 
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by the switch-on function appearing in the noise kernel, 
which eliminates the strong cut-off sensitivities and jolts 
mentioned above when calculating correlation functions 
(the covariance matrix) and its derivatives. 

We conclude this subsection briefly mentioning some 
generally applicable bounds on the growth of the differ- 
ent thermal covariance components. First, we note from 
Eq. (83) that the thermal covariance is positive definite 
as the noise kernel is a positive definite function. We also 
note that the thermal covariance begins with <r T (0) = 
and &t{0) = 0. Given that this matrix is positive def- 
inite, the off-diagonal entries must be smaller than the 
average (arithmetic or geometric) diagonal entries. But 
the off-diagonal cr^p P (t) is proportional to &^ x (t) and we 
have, therefore, the constraint 



< MV /^K P W> (136) 

which is also generally less than the late-time uncertainty 
as both a^ x (t) and <J^(t) begin increasing and then pro- 
ceed to undergo damped oscillations, wherein each cycle 
there is a net increase in uncertainty. This constrains 
the growth of position uncertainty. If the uncertainty 
function takes reasonable values, then the position un- 
certainty must have reasonable growth. 

An analogous constraint can be placed upon the 
growth in momentum uncertainty by considering the pos- 

itive definite matrix $ • v ■ d> which yields 



rPPf 



{t)\ < 2MJaV?(t)(G-v-G)(t). (137) 



So while the growth in position uncertainty is well con- 
strained, growth in momentum is much less constrained. 
Corresponding to this, we show in Sec. IV B 2 that the 
late-time momentum uncertainty has much more sensi- 
tivity to the high frequency modes of the bath. In terms 
of ohmic coupling, the initial linear jolts, &p p ~ A, and 
late-time logarithmic cut-off sensitivity only occurs in the 
momentum uncertainty. The position uncertainty is rel- 
atively well behaved in both respects, having only initial 
logarithmic jolts and no late-time cut-off sensitivity at 
all. The (linear) momentum jolting occurs only for a 
short period of time, At ~ A -1 . The result is a rapid 
momentum dispersion near the initial time, but bounded 
logarithmically. 



Full-Time Diffusion Coefficients for Large 
Cut-off 



Full-time solutions for finite cut-off are completely pos- 
sible given our analytic spectrum, the exact nonlocal 
propagator in Sec. VA, and the contour integrals de- 
tailed in Sec. IV A 2. Such resulting solutions were used 
to plot the early time evolution in Fig. 4, but they are a 
bit cumbersome for publishing. Therefore, for pedagog- 
ical reasons we will restrict ourselves to the high cut-off 



regime in this subsection since substantial additional sim- 
plifications can be employed in that case. For nonlocal 
dissipation it is in general much easier to calculate first 
the thermal covariance than the diffusion coefficients, but 
the situation will be different here. The key point that 
will be exploited in this subsection is that for large cut-off 
the propagator in the ohmic case can be approximated 
by the local one, G R (t), as discussed in Sec. V A. The ad- 
vantage of using the local propagator Gr(£) is that only 
the term involving a single time integral contributes to 
the expression for the diffusion coefficients in Eq. (59). 
On the other hand, if one is only interested in the late- 
time asymptotic values of the diffusion coefficients, one 
can obtain simple analytic results without the need to 
restrict oneself to large values of the cut-off by using the 
results that will be presented in the next subsection. 

The details of the derivation and the complete results 
for the diffusion coefficients at arbitrary times are pro- 
vided in F. Here we simply highlight the main results and 
discuss some of their implications. Both diffusion coeffi- 
cients can be written in the following compact form: 



D xp (t) = D xp (oo) 



(138) 



M 7o { G R (t) + G R (t) (2 7 o-^ ) } DF(f), 



Dp P (t) = D pp (oo 



(139) 



M 7o { G R (t) [lo + j t ) + G K (t) Q 2 } DF(t), 



where D xp (oo) and D pp (oo) are immediately obtained by 
multiplying Eqs. (F3)-(F4) by s and taking the limit 
s — > 0. The general expression for DF(i) is given by 
Eq. (Fll), but a simple result for the zero temperature 
case is provided in Eq. (F14). Essentially, DF(i) decays 
in a manner slightly more complicated than that of expo- 
nential integrals with system, coupling, and temperature 
timescales but such that temperature is the most domi- 
nant. 

It is important to note that the coefficients D xp (t) and 
Dp P (t) both exhibit logarithmic divergences in the limit 
A — > co. This has been pointed out for D xp (t) in Ref. [35], 
where the coefficients of the master equation were cal- 
culated perturbatively to second order in the system- 
environment coupling constants (linear order in 7o ). The 
fact that there is also a logarithmic divergence in D pp (t) 
was not seen in that reference because it is quartic in 
the system-environment coupling constants (quadratic in 
7 o). Moreover, strictly speaking such kinds of perturba- 
tive calculations cannot be employed to study the long 
time behavior since they are only valid for t <C 7o ~ 1 and 
they miss for instance the exponential decay of the second 
and third terms on the right-hand side of Eqs. (F7)-(F8). 

We close this subsection with some remarks about 
the late-time diffusion coefficients in the weak coupling 
regime. Expanding Eqs. (F3)-(F4) perturbatively in 7 o 
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we get 



D. Late-Time Covariance for Finite Cut-off 



D xp (oo) = -70 Re 

7T 







[h(M 









D pp (oo) 



70 f2 coth 



2T 



0( 7o 2 ), 
(140) 

(141) 



In comparison to the weak coupling master equation of 
Caldeira et al. [15], the normal diffusion coefficient is the 
same to lowest order in the coupling, but the anoma- 
lous diffusion coefficient is completely absent there. The 
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FIG. 6. Late time D xp for • high temperature or equivalently 
Caldeira, • HPZ at A = 10 3 fi and A = 10 4 £L 



In Sec. V A the late-time dissipation and renormalized 
frequency coefficients were directly inferred from the non- 
local propagator to be 7* and L\, the result of factoring a 
cubic polynomial in the nonlocal Green function. These 
coefficients are entirely non-perturbative in both coupling 
and cut-off and completely determine the late-time prop- 
agator. The remaining part of the solution pertains to the 
emergence of the thermal covariance, whose late-time dy- 
namics can be described as in Sec. IV B 2, given the late- 
time propagator and the late-time thermal covariance. 
The late-time thermal covariance can also be related to 
the late-time diffusion coefficients through the Lyapunov 
equation, Eq. (68), but the thermal covariance is an eas- 
ier quantity to compute. If interested in the diffusion 
coefficients, one can then obtain them straightforwardly 
using Eq. (68). 

For our spectral density the simplified integrals derived 
in E are contour integrals and can be evaluated via the 
residue theorem after using the rational expansion of the 
hyperbolic cotangent, Eq. (96). The result for the late- 
time, but non-perturbative thermal covariance obtained 
in this way is 



XX 

(Jrp 



m? + 2ik Im[7e] ' 



/run - -/V^f^-i- _ 

T%? = MT + — r-^Im 
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(142) 
(143) 

)■ 

(144) 



largest contribution (in the weak coupling regime) to the 
anomalous diffusion coefficient comes from the cut-off 
and it does not vanish at finite temperature (see Fig. 6). 
This logarithmic sensitivity does not enter into the nor- 
mal diffusion coefficient until second order, but in the 
anomalous diffusion coefficient it is only proportional to 
one power of the coupling constant, which is the or- 
der to which the master equation of Caldeira et al. [15] 
should be valid. In this weak-coupling perturbative ex- 
pansion, both diffusion coefficients are of order 70 plus 
higher-order corrections, but they give contributions of 
different orders to the late-time thermal covariance er^?, 
Sec. IV B 2. Whereas gives contributions of order 1 
because it appears multiplied by a factor I/70, D£2 gives 
contributions of order 70. That is why the correct ther- 
malization in the weak-coupling limit was obtained in 
Ref. [15] despite having completely neglected the anoma- 
lous diffusion coefficient. The origin of the mixed orders 
in 70 appearing on the right-hand side of Eq. (119) can 
be ultimately traced to the fact that H contains terms 
both of order 1 and 70, whose implication for <x^? can be 
straightforwardly seen from the Lyapunov equation (68). 



where we assumed, as before, that fi* = — 7* is real 
and H[z] denotes the harmonic number function defined 
in Sec. Al. If one expands those expressions, and the 
expressions below, under the assumption that L\ is real, 
e.g. using Im[z] = (z — z)/(2z), then one will have the 
more general expressions which will apply even in the 
overdamped regime. 

At high temperature all of the harmonic num- 
ber functions vanish, leaving only the first terms in 
Eqs. (142)-(143), which are proportional to temperature: 



~xx 



MT + O{T ). 



(145) 
(146) 



This corresponds to the high-temperature result of clas- 
sical statistical mechanics. It is interesting that this can 
happen for a finite cut-off and, therefore, outside the 
strict Markovian limit. 

At zero temperature the first terms in Eqs. (142)- 
(143) vanish and all of the harmonic number functions 
can be equivalently evaluated as logarithms, so that 1Z 
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simplified as follows: 

»(^)-^)— ifcMfc) 

+ 0(T). (147) 

This generalizes the results of Unruh and Zurek [16], who 
explored the zero temperature regime in the limit of local 
dissipation. 

Finally, in the weak coupling limit these expressions 
correctly reproduce the free thermal state: 

^-2^ COth (^) + ° (70) ' (W8) 
4 P = ^ coth ^ +O(70) . (149 ) 

One can also see that at weak coupling the uncertainty 
function agrees with the weak coupling approximation 
for moderate values of the cut-off scale, as shown in 
Fig. 7. Had one naively tried to have finite diffusion 
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FIG. 7. Late time AirAp for • high temperature, classi- 
cal statistical mechanics, ■ ■ ■ weak coupling approximation 
| coth ^ ■ HPZ at A = W :i Q. and A = 10 4 SX 

in the limit A — > oo, subtracting by hand the log(A/fi) 
term, one would find a violation of the Heisenberg uncer- 
tainty principle at low temperature and strong coupling 
(see Fig. 8), which renders the theory unphysical. Of 
course this does not happen with the unsubtracted the- 
ory, as seen in Fig. 9. It is, thus, clear that the logarith- 
mic dependence on the ultraviolet cut-off that appears 
in the diffusion is a physically important parameter and 
not something that can be subtracted away. 

While the logarithmic sensitivity appears in both dif- 
fusion coefficients, it is suppressed in the position un- 
certainty by inverse powers of the cut-off. For the mo- 
mentum uncertainty, the logarithmic sensitivity appears 
already to first order in 70 (which is itself quadratic in 




FIG. 8. Late time AxAp for the unphysical, subtracted the- 
ory. 
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FIG. 9. Late time AxAp for the A = 10 3 S1 theory. 



the system-environment coupling constant) and is oth- 
erwise unsuppressed. This behavior had already been 
noticed for Gaussian wave-packets in the Ohmic environ- 
ment [16, 36], and as we have discussed in Sec. IV B 2, 
the position uncertainty will be free of the highest cut-off 
sensitivities for any spectral density. 

Finally, given that our results are nonperturbative, it is 
also interesting to point out what happens in the highly 
nonlocal strong coupling regime mentioned Sec. VA. 
The late-time thermal covariance for this case essen- 
tially corresponds to taking the large f2* limit limit of 
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Eqs. (142)-(143): 



with inverse Laplace transform 
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(150) 
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For this model of strong coupling to the environment, 
and yet finite cut-off, the Brownian particle will become 
strongly localized in position at late time and sufficiently 
low temperatures. And although the particle is localized 
in position, the uncertainty principle is not violated but 
at most minimized in the zero temperature limit. 



VI. SUB-OHMIC AND SUPRA-OHMIC CASES 
A. Sub-ohmic with no Cut-off 

As an example where the nonlocal effects of dissipation 
are important, we will consider one of the most common 
and well-behaved sub-ohmic spectral densities, I(uj) oc 
y/uj, which requires neither a UV nor an IR cut-off in the 
final results (although one still needs to renormalize the 
frequency introducing a logarithmically divergent bare 
countcrtcrm) . Our formulas will take a simpler form if 
we express our spectral density in terms of a quadratic 
coupling constant 7* as follows: 



I(uj) = -M^y/U^U, 



(152) 
(153) 



It is then a straightforward calculation to find the prop- 
agator 



G(s) 



j_ 

M 



s 2 + 2Ty/2^s + O 2 



(154) 



which is amenable to partial fraction decomposition in 
yf s since s is strictly positive. As we have defined our 
nonlinear coupling strength in anticipation of this poly- 
nomial, the roots of the quartic denominator '■ k £ 
{1, 2, 3, 4} can be shown to be the conjugate pairs 



ri,2 



( + y/^ ± lyjw* + 2 7 *) , (155) 

(156) 
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7*3,4 = -7= ( - yfu* ± ~ 2 7* 



After partial fraction decomposition, we may cast our 
propagator in the form 
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(159) 



where erfc(z) is the cumulative error function of the nor- 
mal distribution, defined in A3. There are additional 
terms from the individual Laplace transforms, like t~ x / 2 , 
but they vanish in the sum. Using Eq. (A15) for the 
asymptotic expansion of erfc(z) in order to expand the 
Green function in Eq. (159) at late times, we obtain terms 
of the form 



ze z2 erfc(z) = ^E(- 1 ) 
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(160) 



which we can use to expand the Green function in 
Eq. (159). After grouping all the contributions together, 





FIG. 10. Asymptotic expansion of sub-ohmic ■ propagator 
G(t) into ■ • ■ the local contribution and • the nonlocal con- 
tribution for 7* = 9 . The local contribution is initially more 
significant, but the nonlocal contribution dominates eventu- 
ally. 

we will find exponential terms with characteristic fre- 
quencies / = — T±iy/uj 2 + 27*0;* , which are actually the 
solutions to the characteristic rate equation (114) with 
smallest negative real part. These are the only terms 
that one would have considered if the local propagator 
Gr(<) within the late-time approximation of Sec. IV B 
had been employed. In addition, and more importantly 
are the power-law decay terms which admit no local rep- 
resentation. 

This sub-ohmic model provides a perfect example 
showing when effectively local treatments, such as that 
in Sec. IV B, will fail completely At first the local con- 
tribution will dominate and the master equation coeffi- 
cients will appear to trend towards r as 7* and Or a 
lj* +7*. However, the nonlocal contribution (the power- 
law terms) will eventually dominate the more swiftly de- 
caying local contribution (the exponential terms) and a 
correct treatment of the nonlocal dynamics will be re- 
quired. In fact, as the nonlocal contribution becomes 
comparable to the local contribution, the master equa- 
tion coefficients will become periodically divergent [this 
is related to the fact that det <!?(£) vanishes and changes 
sign at those times.]. The underlying homogeneous evolu- 
tion is well behaved and strictly dissipative (the damping 
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kernel is positive definite), but the localizing perspective 
of the master equation becomes divergently unnatural. 
Any errors, numeric or analytic, can be catastrophic in 
the master equation perspective. In this respect, the sub- 
tleties missed in previous derivations of the master equa- 
tion, as pointed out in Sec. Ill B and which are relevant 
whenever nonlocal effects are important, will likely give 
rise to substantial discrepancies in this case. 

The full-time evolution is rather complicated, but the 
late-time limit is very manageable. For example, from 
Eq. (97) we can express the late-time thermal position 
uncertainty as 

erf? (oo) = 2 / /(w)cothf— )\G(iuj)\ 2 2Vu;dVoj, 
Jo \2ttTJ 

(161) 

where we have used the relation duj = 2^ujdy/uj. The 
integrand is amenable to partial fraction decomposition, 
after a rational expansion of the hyperbolic cotangent 
with Eq. (96), and can therefore be integrated without 
resorting to numerics. Additionally, and in contrast to 
the ohmic case, the integrand is even in yfuj for all tem- 
peratures, including zero, and contour integration tech- 
niques are more generally applicable. 

Strictly speaking we cannot compare exact sub-ohmic 
solutions to those obtained with an incorrect master 
equation since the master equation will yield nonsense, 
but we can compare the exact nonlocal dynamics to those 
obtained by extracting the local dynamics and assuming 
it to be the dominant behavior. Obviously the effectively 
local approximation is incorrect, but it should be good to 
zcroth order in the coupling and one might naively expect 
that it might also behave reasonably for finite coupling 
strength. However, in Fig. 11 we compare the late-time 
uncertainty functions and show there to be sharp dis- 
agreement to the first two orders in the coupling constant 
squared (the slope and the curvature of the curves on the 
plot). 



B. Supra-ohmic with Finite Cut-off 

The conventional wisdom has been to consider supra- 
ohmic spectral densities of the form 



2 , r 



X 



A 



(162) 



where \ '■ [0, oo) — > [1,0) denotes the cut-off regulator. 
Without a cut-off regulator, all supra-ohmic couplings 
have greater than logarithmic high frequency divergence 
in the diffusion and thermal covariance integrals (see E 
for exact integrals in the infinite time limit) . Even when 
regulated, the mere potential for divergence therefore cor- 
responds to cut-off sensitivity from the high frequency 
portion of noise integrals, which is balanced by the ex- 
tra inverse powers of the cut-off in the pre-factor of the 
above spectral density. 
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FIG. 11. Late-time sub-ohmic uncertainty function at zero 
temperature with the • exact nonlocal solution and ■ • ■ ficti- 
tious effectively local solution. In the limit of vanishing dis- 
sipation, one has the minimal uncertainty ground state (zero 
temperature thermal state) in each case. 



Here we will restrict our investigation to the following 
spectral density 



2 (-Y 
'l 



(163) 



because this example is exactly solvable. The corre- 
sponding damping kernel in Laplace space is 



7(s) 



72 



(1 



(164) 



One might be inclined to view this damping kernel as 
a tiny mass renormalization plus even less significant 
higher order terms, but the effect quite different from 
that, as we will see. After factoring the fourth-order 
polynomial, the fully nonlocal propagator can be decom- 
posed by partial fractions into two sets of timescales. Ex- 
panding perturbatively in 72, the first set of timescales 
correspond to the system frequency with weak damping 



7* = 72 



0(7 2 2 ), 
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while the second set of timescales correspond to quickly 
decaying nonlocal contributions associated with the cut- 
off scale: 
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The situation is analogous to ohmic case with a finite 
cut-off except that the nonlocal part of the propagator is 
also oscillating at the rate J7a ~ V72A, for weak coupling 
and high cut-off. 

This form of spectral density was constructed only 
with well-behaved high frequency contributions in mind. 
Nevertheless, as shown in Fig. 12, we find the conven- 





FIG. 12. Late-time supra-ohmic uncertainty function at zero 
temperature for cut-offs between 100f2 and 500Q,. The left 
plot is with a conventional coupling scale, while the right plot 
has decreased the coupling strength by an extra power of the 
cut-off. 

tional form of spectral density to be inadequate. There is 
clearly some cut-off sensitivity in the thermal covariance 
which is remedied by introducing an additional power of 
cut-off suppression. E.g. the conventional form of spec- 
tral density is not well behaved, but the substitution 





72 -> Jl2, 



(169) 



is well behaved. 

An explanation only emerges after a more thorough ex- 
amination of the contour integrals. The high-frequency 
regime, u > A, is already rendered well behaved by 
the conventional cut-off-dependent prefactor. The near- 
resonance regime, u ~ il, which produces the weak cou- 
pling limit, also appears to be well behaved. There is 
only one remaining suspect and it proves to be the cul- 
prit. The previously unaccounted for cut-off sensitivity 
arises here from the nonlocal timescales of the propa- 
gator, i.e. the uj ~ A regime. This is quite surprising 
as unlike sub-ohmic coupling, supra-ohmic coupling does 
yield a well-behaved local representation for its late-time 
dynamics. But residues of the contour integral which 
correspond to the nonlocal timescales reveal the correct 
dominant behavior a pp ~ ^MQ\ = j-M^/t^A, for weak 
coupling and high cut-off. Therefore the conventional, 
linear coupling 72 must be suppressed by an additional 
factor of the cut-off, else the momentum covariance will 
be plagued by a y/X sensitivity. 



VII. GENERALIZATIONS OF THE THEORY 

A. Influence of a Classical Force 

In this section we consider the case of a classical force 
F(t) acting on the quantum oscillator. This is done by 



adding a time-dependent potential —F(t)x to the system 
Lagrangian: 

L s = *M (i 2 - fiV) + F(t)x, (170) 

which gives rise to the following additional source on the 
right-hand side of Eq. (19): 



F(t) 





F(t) 



(171) 



Following our master equation derivation in Sec. IIIB, 
it is easy to see that such a deterministic source in the 
Langevin equation simply adds a driving term to the mas- 
ter equation, which becomes 



d_ 
dt 



W z (z,t) 



(172) 



{ V Z T Uit) z - V Z T Peff(i) + V Z T D(i) V z } W t (z, t) 
where the effective force F c ff(t) is given by 



F eS (t) = F(i) + / dr 
Jo 



*(t-r) F(r) 



(173) 

Note that the last term in Eq. (173) is a consequence of 
having nonlocal dissipation and, as we saw in Sec. IIIB, 
it vanishes for local dissipation. 

Similarly, the method of Sec. II C, based on the solu- 
tions of the Langevin equation, can be straightforwardly 
generalized to this case and one obtains the following 
result for the time evolution of the reduced Wigner func- 
tion: 



W r [t,k] = wJo,$ T (t)k 



e -ik T CTT (t)k e - I k T (z) F (t) 



with a driven mean (z) F (t) given by 



(174) 



(175) 



On the other hand, one can alternatively use the method 
of characteristic curves to solve the master equation, as 
done in Sec. Ill CI. Fourier-transforming Eq. (172), one 
gets an equation analogous to Eq. (66) but with an extra 
term — ik T F c ff(t) on the right-hand side. Following the 
same procedure as in Sec. IIIC 1, one finally obtains the 
same result as in Eq. (174) but with 



(z) F (t) = / dr*(t,r)F cff (r). 
Jo 



(176) 



Eqs. (175) and (176) can be shown to be equivalent as 
follows. First, one rewrites Eq. (173) as 



F off (r) = 



dr 



H(t) 



f dr*(r-T')F(T'). (177) 
Jo 



Next, one substitutes Eq. (177) into Eq. (176) and per- 
forms an integration by parts of the derivative term. 
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Finally, one takes into account that (d/dr) (<& (r)) 
= <fr _1 (r) "H(t), which follows from Eq. (55), and the 
result in Eq. (175) is recovered. Hence, we see that al- 
though <!>(£, t) and <fr(t— r) arc different for nonlocal dis- 
sipation, this is exactly compensated by the contribution 
from the second term on the right-hand side of Eq. (173), 
which does not vanish in that case. 

Note that just as all the temperature dependence ap- 
pears entirely in the second cumulant, or covariance, the 
external force only affects the first cumulant, or mean. 
Eq. (174) shows that the mean, (z)(t), is shifted by 
(z)p{t), which characterizes the response to the driving 
force. In fact, using Eq. (22) one can immediately see 
that it corresponds to shifting (x) and (p) respectively 
by (G * F)(t) and (MG * F)(t), as one would expect. 



B. iV-Oscillator Master Equation 



implied by Eq. (180). After taking the Laplace transform, 
the Langevin equation in position space is then given by 

M (s 2 + 2s7(s) + n 2 ) x(s) = M (s x + x ) + £(s), 

(183) 

which can be solved via matrix inversion to find G(s), 
with which one can construct <&(t) and <TT(t). However, 
closed- form evaluation of G(t) can be rather involved: 
even for local dissipation the two-oscillator problem re- 
quires factoring a fourth-order polynomial. In general, 
the TV-oscillator problem will require factoring a polyno- 
mial of order 2 N for local dissipation and of order 27V + 1 
or more for nonlocal dissipation. We leave more thorough 
discussion to future work, where we will derive block- 
matrix equations for the positon and momentum parts in 
the phase-space representation analogous to those herein. 



Our compact matrix notation allows a number of gen- 
eralizations in a fairly straightforward fashion. As an il- 
lustration we present the generalization of our results for 
the master equation and its solutions to the case of mul- 
tiple system oscillators {x a } (which includes the case of 
a higher dimensional oscillator) with arbitrarily bilinear 
coupling to themselves and to the bath oscillators {yj}- 
We consider the system Lagrangian for N oscillators and 
a generic bilinear term for the system-bath interaction: 

L s = i (x T Mi-xWx) , (178) 

Ut =y T cx = f/V". ( 179 ) 

where we used Einstein's summation convention for re- 
peated indices and the matrix c connects system posi- 
tions (denoted by Greek indices) to bath positions (de- 
noted by Latin indices). The matrices M and MJ1 2 are 
symmetric and positive definite. The eigenvalues of fi 
correspond to the normal-mode frequencies, as can be 
inferred from our Langevin equation (183). 

The effects of the environment for the generalized sit- 
uation described by Eqs. (178)-(179) can be entirely en- 
coded in a simple generalization of the spectral density 
as well as the noise and damping kernels: 

I a p(u,)=Y,8("-"k)** aCkl , (180) 

k, 

v{t,r) = J dwI(w)coth(— ) cos[w(t-<r)] , (181) 

Mj(t,T)= doj^cos[uj(t-T)}. (182) 
Jo u 

In fact, one can directly specify the system-environment 
coupling by giving the spectral density matrix I(w), 
which must be symmetric and positive semi-definite, as 



VIII. DISCUSSION 

Quantum Brownian motion of an oscillator coupled to 
a thermal reservoir of quantum oscillators has been the 
canonical model for the study of open quantum systems 
where one can use it to investigate all the environmen- 
tal effects on an open quantum system it interacts with, 
even of macroscopic scale, such as quantum dissipation, 
diffusion, dccohcrcncc and entanglement. It also pro- 
vides important information on quantum measurement, 
such as noise, fluctuations, correlations, uncertainty re- 
lation and standard quantum limit in mesoscopic sys- 
tems. Many experiments have been carried out for test- 
ing these processes. An exact master equation was re- 
ported some years ago [17] governing the reduced density 
matrix of an open quantum system coupled to a general 
environment of arbitrary spectral density and tempera- 
ture. Subsequently there have also been claims of exact 
solutions [21]. We have found many previous derivations 
to be correct for local dissipation, but containing errors 
or omissions for nonlocal dissipation; in their place we 
have presented the most complete and correct derivation 
of the QBM master equation to date. In this paper we 
report on solutions to this equation for a fairly general set 
of physical conditions and a generalization of the QBM 
master equation to a system with an arbitrary number of 
oscillators. Most of the previous results required one to 
solve integro-differential equations numerically, whereas 
we have reduced everything to quadrature, which can be 
further simplified in many cases using contour integra- 
tion techniques. We expect these results to be useful in 
realistic settings for the analysis of many problems which 
can be described by this model. 

More specifically, we have found a compact expression 
for the general solution of this master equation, show- 
ing that at late times it tends to a Gaussian state en- 
tirely characterized by its asymptotic covariance matrix. 
For odd meromorphic spectral densities, and many oth- 
ers, the result for this late-time covariance matrix can 
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be evaluated as a simple contour integral. As an ex- 
ample we provide explicit exact nonperturbative results 
for an ohmic environment with a finite cut-off which are 
valid for an arbitrarily strong coupling. At sufficiently 
low temperatures and strong coupling this equilibrium 
state becomes highly squeezed and the system becomes 
extremely localized in position space, a phenomenon with 
potentially interesting applications in the realm of meso- 
scopic systems. 

The general solution of the master equation involves 
the matrix propagator of a linear integro-differential 
equation. We have been able to solve these equations 
exactly for several ohmic, sub-ohmic and supra-ohmic 
environments with a finite cut-off and studied the evo- 
lution of the system for finite times. This is achieved us- 
ing Laplace transforms and eventually transforming back 
to time domain. From such exact (and simple) solutions 
for the propagator one gains highly valuable information. 
For instance, one can justify that using the local propaga- 
tor is a valid approximation for the ohmic environment in 
the large cut-off limit. This approximation leads to great 
simplifications and we are then able to provide relatively 
simple analytic expressions for the diffusion coefficients 
of the master equation at all times. Similarly, our ex- 
act solutions for the propagator in specific examples of 
sub-ohmic and super-ohmic environments reveal a dom- 
inant contribution from nonlocal dissipation effects. In 
the first case it is a consequence of long-time correlations, 
due to the low- frequency modes of the environment, that 
become important at late times. In contrast, the source 
of nonlocality in the supra-ohmic case is the UV regula- 
tor function, and it gives rise to a marked cut-off sensi- 
tivity of the momentum covariance which had not been 
noticed so far. On the other hand, it should be pointed 
out that although the results for the exact propagator of 
the integro-differential equation are rather simple, some 
of the general expressions for the solutions of the master 
equation are rather lengthy and have not been reported 
here. They have, nevertheless, been employed to eval- 
uate and plot the exact time evolution of the thermal 
covariance for an ohmic environment with a finite cut-off 
in Sec. VB. 

It is important to discuss the cut-off sensitivity of 
the late-time covariance and diffusion coefficients for an 
ohmic environment in the weak coupling regime. While 
a££ is finite in the infinite cut-off limit, o£2 depends 
logarithmically on A for large A already at order 70 
[Eq. (143)]. This means that it is absolutely necessary 
to consider a finite cut-off. The kind of divergences that 
appear otherwise cannot be dealt with by renormalizing 
the frequency or other bare parameters of the theory. In 
fact, as shown in Sec. VD, subtracting the divergent term 
would lead to inconsistencies (violation of Heisenberg's 
uncertainty principle). Furthermore, from the late-time 
thermal covariance one can immediately obtain the late- 
time diffusion coefficients as well (see the discussion at 
the end of Sec. VC). One finds then that both the normal 
and anomalous diffusion coefficients are logarithmically 



sensitive to large cut-offs. However, while this depen- 
dence appears in D xp [Eq. (140)] at order 70, in D pp it 
only appears at order 7q, and it had been missed in pre- 
vious analytic studies which treated 70 perturbativcly to 
lowest order. 

We would also like to stress the following point. When 
studying an ohmic environment with a finite but large 
cut-off, it can be a good approximation to consider local 
dissipation (infinite cut-off limit for the damping kernel) 
while keeping the cut-off finite in the noise kernel. This 
has already been discussed above and justifies calcula- 
tions like those of Ref. [16] up to corrections suppressed 
by inverse powers of the cut-off. However, the opposite is 
not true: it is essential to keep a finite cut-off in the noise 
kernel to avoid the divergences discussed in the previous 
paragraph. This is precisely the origin of the divergences 
and pathological behavior found in Ref. [21], where a fi- 
nite cut-off was employed in the damping kernel but not 
in the noise kernel. Instead one should use the same 
spectral function everywhere, which means having a fi- 
nite cut-off in both kernels, and everything would be well 
defined then. Note that these divergences would appear 
in the momentum covariance even at asymptotically late 
times, as discussed in the previous paragraph. There is a 
different kind of sensitivity to large values of the cut-off 
that is due to having started with a uncorrelated state 
for the system and the environment. This gives rise to 
a jolt in the normal diffusion coefficient at early times 
of order 1/A with an amplitude proportional to A, as 
well as a logarithmic dependence on the cut-off of a xx 
(and a P p) that decays exponentially with the relaxation 
time-scale 1/r. They would not be present if one had 
started with an appropriately correlated initial sate, and 
then prepared the system in a finite time (not suddenly) . 
Alternatively, this can be implemented by switching on 
the system-environment interaction smoothly in a finite 
time much larger than 1/A, but shorter than the other 
dynamical scales of the system. 

As a further generalization of the QBM master equa- 
tion we have included the influence of external forces. 
This modifies the dynamics by driving the mean position 
and momentum just as with a classical driven system 
(even for nonlocal dissipation). In this model we found 
that the force has no effect upon the width of the wave- 
packet or any cumulant other than the mean. These 
results may be useful for the study of low-temperature 
measurements of driven oscillators, which are relevant 
for experiments with nanomechanical resonators [11, 12]. 
They also play a crucial role in future schemes for the 
detection of gravitational waves with high-intensity laser 



2 A detailed discussion is provided in C 2. There are plenty of 
technical details concerning the initial kick and the effect of the 
switch-on function on the damping kernel, but the real key point 
is the effect of the switch-on function in the noise kernel, when 
evaluating either the diffusion coefficients or the correlation func- 
tions (the covariance matrix). 
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interferometers, where the radiation pressure effects on 
the cavity mirrors are important [37, 38]. 

Finally, we have extended the model of one quantum 
oscillator bilinearly coupled to a thermal reservoir of os- 
cillators to a model of multiple oscillators bilinearly cou- 
pled to themselves and the bath in an arbitrary fash- 
ion. With this generalization, the potential for applica- 
tion [39, 40] becomes almost endless and we leave further 
study to future research [41]. 
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where C( n ) 1S the Riemann zeta function. 

2. Exponential Integral 

The exponential integral is a special function which is 
defined for |arg(z)| < it as 



E x (*) = 




and has a branch cut along |arg(z)| = tt. Its series ex- 
pansion is 

E 1 (z) = - 7E -lnz- V^-z", (A10) 

71=1 

and its asymptotic expansion for \z\ — > oo is given by 

El (z) = ^(l-i + -| + ...Y (All) 



Appendix A: Special Functions 

1. Harmonic Number 

The harmonic number H(n) is a function similar to a 
logarithm, whose definition and main properties are 



1 

H(n) = 1 ' n 6 Z+ (A1) 

k=l 

H(0) = 0, (A2) 
7E = lim [H(n)-log(n))] ; (A3) 

n— >-oo 

where 7e is known as the Euler-Mascheroni constant. Its 
generalization to the complex plane exhibits similar prop- 
erties and is given by 

H(z) = 7E + ^(z + l), zeC, (A4) 

where ip(z) is the digamma function, defined as 

r'(z) 



T(z) 



It satisfies the recurrence relation 



ip(z + 1) = tp(z) + 



(A5) 



(A6) 



and its Taylor expansion around 1 as well as its asymp- 
totic expansion for \z\ — > oo are given respectively by 



k=l 

for \z\ < 1, 

1 1 



tp(z) ~ In z — 



2z I2z 2 
if |arg(z)| < 7T, 



(A7) 
(A8) 



3. Error Function 

The error function is defined as 
2 



erf(z) 



e~ w dw, 



(A12) 



where the path integration is subject to the restriction 
limi, u ,|_ > oo |arg(w) | < ir/4. In addition, the complemen- 
tary error function is defined as 



erfc(z) = — 

The series expansion is 
2 



erf(z) 



e~ w dw = 1 -erf(z). 



\~ L ) z 2n+l 



/?r ^ (2n + l)n! 



and the asymptotic expansion for |z| 
|arg(z)| < 37r/4) is given by 



erfc(z) 



1 



1 3 

2^2 + 4^4 



(A13) 

(A14) 
oo (and 

(A15) 



which along with the fact that erf(z) is odd, is suffi- 
cient to create an accompanying asymptotic expansion 
for |arg(z)| > 37r/4 



Appendix B: Some properties of Laplace Transforms 

Given a real function f(t), defined for all real numbers 
t > 0, its Laplace transform is defined as 

poo 

f(s) = C{f(t)}( S )= / e- st f(t)dt, (Bl) 
Jo- 
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where the one-sided limit from the left for the lower limit 
of integration is chosen so that the transform of the Dirac 
delta function is one, i.e. £{6(t)} = 1. The main prop- 
erties used in the paper are the following. First, the 
Laplace transform of a derivative is given by 



£ 



{/(*)}(*) = «/00-/(o). 



(B2) 



And from this one can easily infer that the Laplace trans- 
form of an integral: 



1 



drf(r)Us) = -f(s) 



(B3) 



Second, multiplying f(t) by an exponential corresponds 
to a translation of the Laplace transform: 



£{e at f(t)}(s) = f( S -a). 



(B4) 



Third, if the inverse Laplace transform of f(s) is /(<) 6(t), 
multiplying /(s) by an exponential corresponds to a 
translation of the inverse Laplace transform: 



£ 



- 1 {e as f(s)}( S )=f(t+a)6(t+a). 



(B5) 



Fourth, the Laplace transform of a Laplace convolution 
is given by the product of the Laplace transforms: 



where 



£{(f*g)(t)}(s) = f(s)g(s), 



(f*g)(t)= / dt'f(t-t')g(t'). 
Jo 



(B6) 



(B7) 



Fifth, the initial value theorem relates the initial value 
of a function f(t) and the infinite limit of its Laplace 
transform as follows: 



/(0 H 



lim sf(s) 



(B8) 



Bromwich's integral illustrates the close relationship 
between the Laplace transform and the Fourier transform 
through analytic continuation. However, even if all the 
singularities of /(s) lie on the Re(s) < half of the com- 
plex plane, the relation is not direct because the Laplace 
transform involves an integral with domain [0, oo) rather 
than (—00,00). The precise relationship can be under- 
stood as follows. Consider a real function f(t) defined for 
all real values of t and whose Fourier transform is f{uj). It 
is useful to define the following additional Fourier trans- 
forms: 



+00 



dtt 



t f(t)9(±t). 



(Bll) 



such that J(uj) = + /-(w) and which satisfy the 

property f±(—oj) = (/±(w)) since f(t) is real. Assum- 
ing that the Laplace transform f(s) has no singularities 
for Re(s) > 0, it can be related by analytic continuation 
to f+(u): 



/+M = hm/(e-Hw). 

e— )-0 



(B12) 



If f(t) is an even function, one has = /+(— uj), and 

using Eq. (B12) one can then write 



f(u) = f+(u) + = lim /(e+iw) + /(e-tw) 

c— >0 



(B13) 



Similarly, if f(t) is an odd function, one has /-(w) 
—/_!_(— u), which implies 



/M = /+M - = lim f(e+iu) - f(e-iw) 

e— >0 



(B14) 



Appendix C: System-Environment Interaction and 
Renor malizat ion 



Sixth, the final value theorem relates the infinite limit 
of a function f(t) and the initial value of its Laplace 
transform as follows: 



/(oo) = lim sf(s), 



(B9) 



provided that all the poles of f(s) are on the Re(s) < 
half of the s complex plane. 

Seventh, the inverse Laplace transform of f(s) can be 
calculated using Bromwich's integral, which involves an 
analytic continuation of f(s) in the complex plane: 

-1 pa+too 

f(t)=£- 1 {f( S )}(s) = — e st f(s)ds, (BIO) 



1. Frequency Renormalization and the Damping 
Kernel 



In this subsection we discuss the relationship between 
the Lagrangians in Eqs. (1) and (2). The question of 
frequency renormalization plays a central role in this dis- 
cussion and it will be analyzed by rewriting the equa- 
tion of motion in terms of the damping kernel, given by 
Eq. (8), rather than the dissipation kernel. Here we will 
take 9 s (t) = 1 and leave the examination of effects due to 
a non-vanishing switch-on time for the next subsection. 

If one starts with the Lagrangian in Eq. (1), the homo- 
geneous part of the Langevin integro-differential equation 
analogous to Eq. (7) is then 



where a is a real number chosen so that the integration 
path lies within the region of convergence of /(s), i.e., 
a > Re(sj) for every singularity Sj of f(s). 



(L-x)(t) = Mx(t) + MD,i are x(t) +2 / dr fi(t-T)x(r). 

(CI) 
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Using the expression for the dissipation kernel in terms 
of the damping kernel, /i(t— r) — M(d/dt)j(t—r), and 
integrating by parts, one obtains 



(L-x)(t) =Mx(t) + 2M dT^(t~r)x(T) 



(C2) 



+ M(nl aie -Sn 2 ) x(t) + 2M 7 (t) x(0), 

where the time-independent frequency renormalization 
<5f2 2 is given by 



sn 2 



2 7 (0) 



2 

M 



dco 



/(«) 



(C3) 



By choosing a bare frequency with an appropriate coun- 
terterm in order to cancel the frequency renormalization, 

^barc = n2 + <^ 2 > one would be finall y left with E q- ( 9 ) 

and an effective frequency for the system oscillator Q. 
This is the same equation that one obtains if one starts 
with the Lagrangian in Eq. (2) and takes 9 s (t) — 1, which 
can be easily understood as follows. Recalling the defi- 
nition of the spectral function I(lo) in terms of the cou- 
pling constants c„ , it is immediate to see that the square 
of the last term on the right-hand side of Eq. (2) gives 
— (l/2)M(50 2 x 2 . Therefore, the Lagrangians in Eqs. (1) 
and (2) are equivalent provided that one makes the choice 
mentioned above for f^ 2 ,™- 



2. Initial-Time Divergences, Coupling Switch-on 
and Initial-State Distortion 

a. Initial-Time Divergences and Coupling Switch-on 

The derivation of the HPZ master equations relies 
upon the key assumption that the system and environ- 
ment are initially uncorrelated. For an ohmic environ- 
ment, this gives rise to an initial "jolt" in the normal 
diffusion coefficient of the master equation with a char- 
acteristic time-scale of order A -1 and an amplitude pro- 
portional to A. Similarly, the frequency ^^(i) in the 
master equation starts with a large value of the order of 
A and decreases to moderate values in a time of order 
A- 1 . 

The physical origin of the jolts in the coefficients of 
the master equation as well as other initial time diver- 
gences, such as the divergent contributions to correlation 
functions of system observables that are due to diver- 
gent boundary terms at the initial time (see Appendix D 
in Ref. [ ]), can be understood as follows. In general 
when a system couples to an environment with an infi- 
nite number of modes, well-behaved states exhibit cor- 
relations with arbitrarily high-frequency modes. In con- 
trast, states that are uncorrelated for sufficiently high 
frequencies (such as completely factorizable states) are 
pathological. For instance, in the limit of infinite cut-off 
they have infinite energy (even with an origin of ener- 
gies such that the ground state of the whole interact- 
ing system has vanishing energy) and their Hilbert space 



is unitarily inequivalent to the space of physical states, 
spanned by the basis of energy eigenvectors of the whole 
system Hamiltonian including the system-environment 
interaction. (Of course for a finite UV cut-off there are 
no divergences or unitary inequivalence, but the poten- 
tially divergent terms are very sensitive to changes in 
the value of the cut-off.) Physically acceptable initial 
states that correspond to the thermal equilibrium state 
for the whole system can be obtained using Euclidean 
path integrals [43]. However, the instantaneous prepa- 
ration functions employed in Ref. [43] to produce other 
states in addition to the thermal equilibrium state still 
give rise to initial divergences, as explained in Ref. [44]. 
In order to obtain finite results, one needs to prepare the 
new initial state within a non- vanishing time [34], which 
corresponds to a physically more realistic situation. The 
alternative approach that we follow here is to switch on 
the system-environment interaction smoothly within a 
time t s much longer than A -1 but shorter than any other 
relevant time-scale of the problem. In this way the fac- 
torized initial state, which is perfectly acceptable in the 
uncoupled case, becomes adequately correlated with the 
arbitrarily high-frequency modes in a regular fashion. 

When adding the short time switch-on function to the 
system-environment coupling to turn on the interaction 
gradually, as in Eq. (2), the initial jolt is no longer present 
in the results for the diffusion coefficients, which behave 
smoothly during the switch-on time. Furthermore, for 
times much longer than t s the contribution to Eq. (59) 
from the switch-on period is negligible and one can sim- 
ply use that equation without including any switch-on 
function. This point is implicitly exploited throughout 
the paper: unless explicitly stated, our calculations of the 
diffusion coefficients do not take into account the switch- 
on functions and the results for those coefficients should 
only be regarded as valid for times sufficiently larger 
than i s , while their values during that period should be 
smoothly interpolated so that they vanish at the initial 
time. 

Either the quick transition from the bare frequency to 
the renormalized one (in the absence of a smooth switch- 
on function) or switching on the interaction in a finite 
time t s can have a non- negligible effect on the homoge- 
nous solutions of the Langevin equation even for times 
much larger than A -1 or t s . Fortunately, as we will show 
in the remaining subsections, the effect can be entirely 
accounted for by a finite shift of the initial momentum 
and the corresponding transformation of the initial state. 



b. Initial Kick (finite cut-off, vanishing switch-on time) 

We start by considering the case in which there is no 
switch-on time and analyze the effect of the slip term, 
which corresponds to the last term on the right-hand 
side of Eq. (9) for the Langevin operator. It can be in- 
terpreted as a transient driving force 



FJt) = 2 7 (i)x , 



(C4) 
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whose contribution to the solution is simply adding a 
term G(t) * Ey(t). The infinite cut-off limit can be ana- 
lyzed using distributions in the time domain or working 
with Laplace transforms. As derived in Eqs. (13)-(14), 
the solution of Eq. (9) (including the slip term) is given 
in Laplace space by 

x(s) = M (sx + x Q ) G(s) + G(s)i(s), (C5) 

whereas one can easily infer that the solution without the 
slip term would be 

y(s) = M ( S y + yo + 2j(s) y ) G(s) + G(s)i(s). (C6) 

In the limit of local dissipation (large cut-off limit) 7(5) = 
70 and one can see that the effect of the slip term is an 
initial kick to the homogeneous solutions, whose values 
before and after the kick can be related via 

yo = x , (C7) 
Vo = i - 2j x . (C8) 

This induces a distortion of the reduced Wigner func- 
tion associated with the transformation xq — > xq — 270X0 
which occurs within the cut-off timescale. The effect of 
such an initial kick can be entirely absorbed in a redefini- 
tion of the initial state, as will be discussed in Sec. C 2 d. 

c. Initial Kick (large cut-off, non-vanishing switch-on time) 

Next, we consider the case with a non- vanishing 
switch-on time t s and smooth switch-on function such 
that 6 S (Q) = and 6 s (t) = 1 for t > t s . Integrating the 
dissipation kernel by parts in Eq. (7), the homogeneous 
part of the Langevin equation becomes 

x(t)+2f dTj F (t,T)x(r) + n 2 x(t) = 
Jo 

-20 B (i) [ dT>y(t-T)6 s (T)x(T), (C9) 
Jo 

in terms of the positive-definite kernel 

lP (t,r)=j(t-r)e a (t)e s (r), (CIO) 

where we have not made any approximations concern- 
ing the timescales of the dissipation kernel and switch-on 
function yet, but have expressed our result in terms of 
the damping kernel defined in Eq. (8) and taken into ac- 
count Eq. (C3). Either in the limit of local dissipation 
or vanishing switch-on time, the term on the right-hand 
side of Eq. (C9) gives rise to a slip term analogous to that 
found in the previous subsection. This is because for x(t) 
evolving slowly compared to A -1 or t s , one has a a con- 
volution of the distributions j(t) and 6 a (t), which is also 
a distribution localized near the initial time. In partic- 
ular, it is immediate to see that the results of Sec. C 2 b 
are recovered in the limit of vanishing t s [see the remark 
below about 9 s (t) in that limit]. 



If we take the high cut-off limit 7(i— t) = joS(t— r), 
which should be a good approximation for A>f~ 1 , the 
right-hand side of Eq. (C9) takes a simple form and we 
are left with 

x(t) + 2 7o 9 2 s {t) x(t) + n 2 x(t) = -70 S s (t) x(t) (Cll) 
5 s (t) = ~el{t), (C12) 

where 8 s (t) is a representation of the delta function in 
the limit of vanishing switch-on time; however, its sup- 
port is entirely contained in the t > interval, so that 
J °° S B (t)dt — 1. [We also took the local dissipation limit 
on the left-hand side of Eq. (Cll) for simplicity, but we 
needn't have done so: that is a slowly varying term which 
does not play an important role here.] For a very rapid 
switch-on function we have 6 s (t)x(t) « 6 s (t) x(0) and 
this term produces an initial kick, xq — > xq — 70^0 ; anal- 
ogous to that described in the previous subsection. This 
kick of the homogeneous solutions will produce a distor- 
tion of reduced Wigner function which occurs within the 
switch-on timescale. For times much larger than t s , this 
effect can also be entirely absorbed into a redefinition of 
the initial state, as described in the next subsection. 



d. Initial-State Distortion 

In Sec. C 2 b we calculated that in a particular limit 
of £1 <C A <C tjT 1 one obtains a kick to the initial state 
of Xq — > Xq — 270X0 which occurs within the slower cut- 
off timescale. Whereas in Sec. C 2 c we calculated that 
in a particular limit of -C t^ 1 -C A one obtains a 
kick to the initial state of xq xq — 70X0 which occurs 
within the slower switch-on timescale. From the exact 
relation in Eq. (C9), if one tries to enforce both high 
cut-off and short switch-on time then there will be a kick 
Xo — > Xq — cjqXo which occurs in the slower of the cut-off 
and switch-on timescales. And if the stationary damp- 
ing kernel j(t) and switch-on function's derivative 6 s (t) 
are suitably well-behaved distributions, then this kick is 
bounded so that < c < 2. 

From these results one might be tempted to con- 
sider modifying the Lagrangian by introducing a suitable 
time-dependent frequency renormalization counterterm 
^kickM = ~ c laS(t)- However, even though an appro- 
priate choice of time-dependent counterterm could com- 
pensate and effectively remove the effect of the initial 
kick in either case, a truly finite cut-off is still necessary 
to have a finite thermal covariance, and the switch-on 
function for the system-environment interaction is still 
essential to avoid the highly cut-off sensitive initial jolt 
in the normal diffusion coefficient and other irregulari- 
ties associated with an uncorrelated initial state [the key 
point in these cases is the dependence on the switch-on 
function of the noise kernel shown in Eq. (5)]. 

Moreover, the effect of any such kick can easily be ac- 
counted for by simply distinguishing between the "bare" 
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initial state before the kick and the "renormalized" state 
immediately after the kick. Following the approach in 
Ref. [19] one can easily see that this initial kick trans- 
lates into a distortion of the Wigner distribution from 
the bare initial state to a shifted one 

Wbaxe(x,P) -> W rcn (x,p) = W hal - c (x,p-cMj x)- (C13) 

This phase-space transformation has a Jacobian matrix 
K with determinant equal to one: 



K 



1 

-cMjo 1 



dctK = 1. 



(C14) 



Therefore, it is simple to calculate renormalized expecta- 
tion values in terms of bare expectation values and vice 
versa: 



(A(x,p))™ = dxdpA(x,p)W™(x,p), (C15) 

bare J J bare 

(A(x,p)) Ten = {A(x,p+cM ln x)) h&IC . (C16) 

We can immediately see that the normalization, linear 
entropy (see Sec. IV A 3) and state overlap are all un- 
changed by the kick. We can also check explicitly that 
the Heisenberg uncertainty relation is preserved as fol- 
lows. First, we start with the covariance matrix for x 
and p corresponding to the Wigner distribution 



' xx @xp 



(T 



with o . 



px &pp 
&px 



(C17) 



(xp) 



ren and Cpp 



pp — (pp)ron, and which transforms in the following way 
er linear phase-space transformations: 

(C18) 



ro- 
under linear phase-space transformations 

er -> KcrK T . 

Hence, from Eq. (C14) we have 

det cr barc = det cr ron . 
Finally, one takes into account that 

h 2 



(det er) > 



4 ' 



(C19) 



(C20) 



corresponds to the formulation in terms of the Wigner 
function of the generalized Heisenberg uncertainty rela- 
tion due to Schrodinger [45, 46]: 



(Ax 2 ) (A P 2 ) 



1 \ ti 2 

-{Ax,A P }\ >-, (C21) 



where {A, B} = AB + BA. 

Furthermore, by switching to the density matrix for- 
malism, we can see that pure states are mapped to pure 
states and positivity is preserved. It is a straightforward 
calculation to show that 



/2bare(*£; V) ^ Prcn 



(C22) 



Pxen(x,y) 



- m to -2 



cMtq 2 



Therefore, if we start in a pure state, which acts as a 
projection operator 



Pba 



(C24) 



then it is fairly easy to see that this will hold for the dis- 
torted state. Additionally, given the positivity condition 



|pbarc#) > 0, 



(C25) 



for all vectors then it is easy to see that the distorted 
state will also fulfill this condition by simply considering 
the vectors e lcMl ° x / 2 i(;(x) in position representation. 

In summary, the new Wigner function that results from 
the transformation defined by Eq. (C13) always corre- 
sponds to a physical density matrix since the transforma- 
tion preserves the normalization and the real-valuedness 
of the Wigner function (implying the normalization and 
hermiticity of the density matrix) as well as the pos- 
itivity of the associated density matrix. Therefore, if 
one is interested in analyzing the evolution of a certain 
state of the system better correlated with the environ- 
ment, one can simply take such a state as W rcrL (x,p) and 
study its evolution for t ^> max[t s ,A -1 ] by considering 
the Langevin equation without the term that gives rise 
to the initial kick. However, given any W ren (x,p) it is 
always possible to follow in detail the evolution during 
the switch-on time by inverting Eq. (C13) to obtain the 
corresponding initial Wigner function before the interac- 
tion was switched on and using the full Langevin equa- 
tion with the contribution from the right-hand side of 
Eq. (C9) included. In general this approach can be re- 
garded simply as a formal procedure to generate a bet- 
ter correlated initial state, but the explicit construction 
involving unitary evolution for the whole system at all 
times guarantees that the result is well defined (i.e. the 
exact solutions of the master equation obtained in this 
way preserve the positivity of the density matrix). 3 



Appendix D: Peculiarities of Propagators and Green 
Functions Associated with Integro-differential 
Equations 

In this Appendix we discuss a subtle mathematical 
point which, to the best of our knowledge, has been 
missed in the existing literature on master equations of 



x P^{x,y)e- l - v . (C23) 



3 Using this approach the system-environment correlations at high 
frequencies will be the same as those of other properly corre- 
lated states (such as the global equilibrium states considered in 
Ref. [43] or states prepared from those in a finite time). However, 
in general the correlations for low frequencies will differ and the 
states of the whole system plus environment will not be equiv- 
alent even if their reduced Wigner functions are the same. In 
particular this implies that even if the reduced Wigner functions 
of the two states coincide at some given time, they will in gen- 
eral evolve differently (until thermal equilibrium for the whole 
system is reached). 
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QBM models. It has to do with properties of Green func- 
tions which are satisfied for ordinary differential equa- 
tions but not for integro-differential equations. Thus, it 
becomes particularly relevant whenever the nonlocal as- 
pects of the dissipation kernel cannot be neglected. 
Consider an integro-differential equation of the form 



z(t) 



with the kernel H(t- 
can be written as 



rfrH((-T)z(r) = £(t), 



(Dl) 



t) given by Eq. (20). Its solutions 



z(t) = *(i)z + (**£)(*), 



(D2) 



where Zq specifies the initial conditions and the matrix 
propagator <&(i) is given by Eq. (24). As far as the homo- 
geneous solutions are concerned, the values of a solution 
at two different times r and t are related by the transi- 
tion matrix $ _1 (t). On the other hand, for some 
given initial conditions the inhomogeneous solutions are 
obtained by integrating the source with the retarded ma- 
trix propagator <& rot (t— r) = $(t— t) 6(t— r), as shown in 
Eq. (D2). 

In the case of a linear differential equation (i.e. for a 
local damping kernel), the retarded matrix propagator 
and the transition matrix are related in a simple way: 
&ret(t— t) = <&(<) 4> _1 (r) 0(t— r). This can be seen by re- 
alizing that $>(t) 0(t—r) satisfies the differential equation 
except for a delta function that results from differentiat- 
ing the theta function, and that the two expressions are 
equal to the identity matrix at t = r. In contrast, for an 
integro-differential equation (a nonlocal damping kernel) 
&(t) $ _1 (t) 9(t—r) no longer corresponds to the retarded 
matrix propagator because <fr(i) 9{t— r) does not satisfy 
the integro-differential equation, which can be seen (for 
t > 0) as follows: 



*(*) 



dr'H(f-r')*(r') + 



dr'H^-r')*^'), 
(D3) 



where the right-hand side equates to 

-J dr'H(t-T')&{T') = -J dr' ' H(i-r') *(t') 0(t'-t). 

(D4) 

The discrepancy is due to a term of the form J^dr' H(i— 
t') 3>(t') (with t > r), which vanishes in the case of non- 
local damping kernel and hence a nonlocal kernel H(i— r'), 
but does not vanish in the nonlocal case. On the other 
hand, 4> rot (t—r) does satisfy the integro-differential equa- 
tion with a delta source, as it should. This point, which 
can be alternatively seen in Laplace space fairly easily, 
follows from the fact that <fr(i) is a solution of the integro- 
differential equation by construction, together with the 
translational invariance of this kind of solutions [i.e. if 



<!?(£) is a solution, <!>(£ — r) is also a solution 4 ]. Such 
a translational invariance follows quite straightforwardly 
from the causal and translationally-invariant nature of 
the kernel H(i — t') as well as the matrix propagator's 
support only for non-negative values of its argument: 



*ret(t-T) 



-/ dT'H(t-T-T')<f>(T')+I5(t-T) 

Jo 

-J dT"H(t-T")$(T"-T)+I5(t-T) 

-f dr" H(t-r") #(r"-r) +IS(t-r), 
Jo 



(D5) 

where I is the identity matrix and we used the fact that 
3?(t') = for t' < in the last equality. 

From the previous discussion it immediately follows 
(taking t > r) that, contrary to the local case, the matrix 
propagator does not factorize in the nonlocal case, i.e. 



*(t-r) ^ *(t)* _1 (r). 



(D6) 



This lack of factorizability also implies that the Green 
function or, equivalently, the matrix propagator <&f(i,r) 
for the integro-differential equation when the boundary 
conditions are specified at some final time [and given by 
Eq. (28)] is no longer an advanced propagator, i.e. it is 
no longer true that $f(i,r) = for t > r. This can be 
proved by contradiction. If one considers r > t > t' in 
Eq. (28) and assumes that 3?f(r, r') = 0, one is left with 







-*(r,f)*(t— 7-') + *(t-t')- 



(D7) 



Taking the limit t' t~ of Eq. (28) and taking into 
account that lim M _ i .o+ <&(u) = I, one finally obtains 
$(r — t) = $(r)$ _1 (i), which is in contradiction with 
Eq. (D6). Therefore, the assumption $f(r, r') = for 
t > t 1 cannot be true in the nonlocal case. 

These facts or closely related ones have been missed 
in the existing literature on master equations for QBM 
models. As a consequence, the existing results for the 
coefficients of the master equation are mathematically 
incorrect unless strictly local dissipation is considered, 
and can give rise to significant discrepancies whenever 
nonlocal effects are important. We close this appendix 
by briefly describing how this affects the different exist- 
ing approaches to deriving the exact master equation for 
QBM models. One class of derivations [18-20] involve 
an intermediate step where the solution of an integro- 
differential equation like Eq. (Dl) with specified bound- 
ary conditions (position and velocity) at a final time is 
needed. The previous discussion directly applies to this 
class of derivations and the main consequences are that 



4 Note that if one uses a convention according to which $ (t ) =0 
for t < 0, then the notation 0(t) is redundant. 
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the Green functions appearing there are not advanced 
and the explicit expressions which were provided, based 
on the assumption that those Green functions were ad- 
vanced, are incorrect. Nevertheless, the results in those 
references can be easily corrected by removing the quali- 
fication of "advanced" propagator and discarding the ex- 
plicit expressions for that Green function. The results 
would then become equivalent to the general result that 
we have obtained in Sec. IIIB, although one would need 
to find a way to construct the Green function explicitly. 
We provide such an explicit construction of the corre- 
sponding matrix propagator $f(r, r') in Eq. (28), where 
it is expressed in terms of known quantities, namely, <f*(t) 
as given by Eq. (24). Note, by the way, that if one 
had truly advanced propagators, one could show that 
the terms involving triple time integrals in the results 
for the diffusion coefficients [such as Eqs. (B.17)-(B.18) 
in Ref. [19]] actually vanish. In fact, these terms cor- 
respond to the last term on the right-hand side of our 
Eq. (51), which only vanishes for local dissipation, as can 
be seen from Eqs. (59), (55) and the discussion above. 

A second class of derivations, including HPZ's original 
derivation for arbitrary temperature and spectral func- 
tion, relies on the use of Green functions for the same 
integro-diffcrcntial equation, but associated with mixed 
boundary conditions which correspond to specifying the 
position at the initial and final times. Explicit expres- 
sions are provided for those Green functions G(t, s) in 
terms of homogeneous solutions ui(t) and U2(r) which 
vanish at the final and initial times respectively. Unfor- 
tunately, although those expressions are standard results 
for ordinary differential equations, they are not valid for 
nonlocal integro-differential equations. This is because 
they involve the sum of two terms, each one of them 
being a certain solution of the homogeneous integro- 
differential equation times 9(t—s) and Q(s—t) respectively 
[see Eq. (2.34) in Ref. [ ]]. However, for similar reasons 
to those given above and illustrated by Eq. (D3), when 
multiplied by the theta functions those solutions cease to 
satisfy the integro-differential equation. 

Finally, a third class of derivations [ ] are based 
on showing that the solutions of the Langevin equa- 
tion can be equivalently understood as solutions of a lo- 
cal ordinary differential equation rather than an integro- 
differential one. This is true for the homogenous solutions 
of the Langevin equation and corresponds to the equiva- 
lence (after inverting and transposing) between the ma- 
trix propagator $(i) associated with the Langevin equa- 
tion and the matrix propagator *frfc(£) associated with 
the ordinary differential equation (72), which we found in 
Sec. IIIC 1. However, such an equivalence is not true for 
inhomogenous solutions of the nonlocal Langevin equa- 
tion. One way of seeing this is by realizing that since 
Eq. (72) is an ordinary differential equation, its retarded 
matrix propagator does factorizc. But if the inhomoge- 
neous solutions of the local equation constructed with 
that propagator were also solutions of the inhomoge- 
neous Langevin equation, it would imply that the re- 



tarded propagator associated with the latter also fac- 
torizes, which is not true for nonlocal dissipation, as 
we showed above. ° In particular, the derivation of 
Eq. (2.18) in Ref. [21] is valid if one takes a vanishing 
inhomogeneous source F(t). Nevertheless, when deriv- 
ing Eq. (2.18) for a non-vanishing source, the authors 
implicitly assumed that if the homogenous solutions of 
the Langevin equations satisfy a local differential equa- 
tion, the inhomogeneous solutions of the Langevin equa- 
tion should also satisfy the inhomogeneous version of the 
same local equation. As we have explained, it turns out 
that this is only true for local dissipation. Not surpris- 
ingly, making use of Eq. (2.18) the authors derive a mas- 
ter equation with diffusion coefficients lacking the terms 
with triple time integrals mentioned above, which in re- 
ality should only vanish for strictly local dissipation. 



Appendix E: Derivation of the Late-Time Thermal 
Covariance 



Here we present the derivation of the general single- 
integral representation of the late-time thermal covari- 
ance. For the sake of brevity we will work out the ex- 
plicit case of the late-time position uncertainty. The 
late-time momentum uncertainty is analogous and the 
cross-correlation vanishes at late times, as implied by 
a* p (t) = (M/2)af : (t) if af(t) tends to a constant 
asymptotic value. 

We start with the full-time, exact expression 



of (t)= J dw/(w)coth(^) (El) 

X dTl dT 2 G(T 1 )cOs[uj(T 1 -T2)}G(T2), 

Jo Jo 

where we have made the simple change of variables t[ — 
t — Ti for i = 1,2. Introducing the additional change of 
variables r = T\ + r 2 , the result can be rewritten as 



oo 

off B (t) = J du/Hcoth^J (E2) 

ft rT 2 +t 

x dr 2 dfG(f-T 2 ) cos[w(f-2 r 2 )] G(r 2 ). 

JO Jt 2 

The double time integration can then be split into two 
parts: 

dr 2 df= I dr 2 d? + dr 2 df. (E3) 
Jo Jt 2 Jo Jt 2 Jo Jt 



(w)coth^j 



To use this argument directly one should consider the equation 
satisfied by [4>^(t)] _1 rather than Eq. (72), which is satisfied by 
#fc(t). That equation can be easily obtained by transposing and 
taking the matrix inverse of Eq. (72) applied to #^(t), and it is 
still a local linear differential equation. 
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At sufficiently late times the contribution form the second 
integration domain can be neglected and we can approx- 
imate the whole integral as follows: 



dr 2 df « / dr 2 df = df dr 2 , (E4) 
Jo Jt 2 Jo Jt 2 Jo Jo 

The next step is to express the cosine in complex form 
with exponential functions. Once that is done, it is 
not difficult to manipulate the result into the form of 
a Laplace convolution: 



*r(t)« fdw/M «>th(^) 



(E5) 



x / drRe{[e- MT G(r)] * [e+^ r G(r)] } , 



where we renamed r as r. Using the property of fre- 
quency shifting in the Laplace domain, i.e. £{e A * f(t)} = 
f(s— A), we obtain 

&j?(s) « / du}I(ui) coth( — J - G(s+zu;) G(s— jcj). 



(E6) 



Application of the final value theorem, as given by 
Eq. (B9), then immediately reveals the exact late-time 



Appendix F: Moderate-Time Diffusion for Ohmic 
Case with Large Cut-off 



In this appendix we calculate the diffusion coefficients 
for the ohmic case using the local propagator Gr(£) in- 
stead of the exact one, which is a valid approximation in 
the high cut-off regime, as discussed in Sec. V A. The big 
advantage of using Gr(*) is that only the first term on the 
right-hand side of Eq. (59), which involves a single time 
integral, will give a non-vanishing contribution. Further- 
more, the Laplace transforms of the corresponding equa- 
tions for the diffusion coefficients exhibit a rather simple 
form if one takes the following steps. First, one writes 
the cosine of the noise kernel in exponential form; next, 
manipulates the time integral until one has a Laplace con- 
volution; and then uses frequency shifting in the Laplace 
domain, i.e. e xt f(t) — > f(s — X). After some algebraic 
manipulations one finally gets 



- ' J d^I{uj) coth^^^Re G R (s+iu;) 



D xp (s) 



(Fl) 



1 



D pp (s) = +- / dw/(w)coth( — Re Gr(s-Nw) 



2T 



(F2) 



covariance 



c4 x (oo) = j duI(uj)coth(^G{+iLj)G(-iLj). (E7) 

Proceeding in a completely analogous way, one can ob- 
tain the result for the momentum covariance and the 
cross correlation. For the cross correlation, the time 
derivative of one of the propagators gives an extra fac- 
tor (s + iui) in the expression in Laplace space. When 
taking the real part, as in Eq. (E6), one is left only with 
s, which cancels out the factor 1/s in Eq. (E7). Appli- 
cation of the final value theorem, as given by Eq. (B9), 
gives then a vanishing result for the asymptotic value of 
the cross-correlation: a^F(oo) = 0. As for the momen- 
tum covariance, the two time derivatives, one for each 
propagator, give an extra factor (s 2 + u> 2 ) in the expres- 
sion in Laplace space. When taking the real part and 
applying the final value theorem, one is left with 



<r£ p (oo) = M' 



duj I (oj)coth(^-^j oj 2 G(+iuj)G(-illi). 

(E8) 



Our late-time Green function (112) is rational in 
the Laplace domain [with late-time coefficients given 
by Eq. (132)]. Moreover, the spectral density I(u>) in 
Eq. (123) is meromorphic with a finite number of poles. 
Together with the rational expansion of the hyperbolic 
cotangent in Eq. (96), this implies that the frequency in- 
tegrals over to in the above diffusion coefficients become 
sums over k of trivial contour integrals in the Laplace 
domain. Still in the Laplace domain, these sums can be 
identified as harmonic number functions (or, equivalcntly, 
digamma functions): 6 



D xp (s) 
Dpp(s) 



_2to_T 

As 
270T 
s 



7o 



ImpZs 



s \ 
A 



To 



Im 



7o- 



(F3) 
, (F4) 



Taking into account Eqs. (E7)-(E8) and the vanishing 
value of the asymptotic cross correlation, the asymptotic 
value of the thermal covariance matrix can be written as 

<T T (oo) = Sy / dw/(w)coth(— J M+^)ioi}& (-«•>), 

(E9) 



where Sy denotes matrix symmetrization. 



6 Many of the expressions derived throughout this paper assume 
underdamping, i.e. 70 < fl with Q = ^/f2 2 — 7q- They can be 
used for the overdamping regime by making the following analyt- 
ical continuation: f2 



17 with 7 = W 7q — f2 2 real. Therefore, 
Eqs. (F3)-(F4) can be applied to the overdamping case if the Im 
and Re terms are first expanded assuming that tt is real, e.g. 
using Im[z] = (z — z)/(2i), and then the analytical continuation 
Q — > 17 is made. 
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in terms of the dimensionless quantities T s and T s denned 

-l 



Is 




-H 



(F5) 



(F6) 



and where 7 S = 70 + s. Note that by making use of 
the final value theorem in Eq. (B9), we only need to 
discard the overall 1/s factor and replace 7 S with 70 in 
Eqs. (F3)-(F4) to obtain the late-time asymptotic values 
D xp (oo) and D pp (po). The H(z) functions are the har- 
monic number function discussed in Al. These terms 
make up, among other things, the well known log(A/f2) 
divergence. They behave asymptotically like logarithms 
but with H(0) = 0, making both their high and zero 
temperature limits trivial. At high temperature, all of 
the harmonic number functions vanish, leaving only the 
second terms which are proportional to the temperature. 
At zero temperature, all of the harmonic number func- 
tions can be equivalently evaluated as logarithms. 

The diffusion coefficients can be expressed in the time 
domain as their asymptotic values plus damped oscillat- 
ing differential operators acting on the same decay func- 
tion DF(i) (although the sums over k cannot in general 
be identified with any simply behaved special functions) : 



D xp {t) = D xp (oo) 



(F7) 



M 7o { G R (t) + G R {t) (2 7o -- ] } DF(i), 



D pp (t) = D pp (oo) 



(F8) 



M l0 { G R (t) ( 7o + ^ ) + G R (t) n 2 } DF(t), 



with the thermal decay function 
cot (±) 



DF(i) = - 



(i + i) 2 + (S 



;TS(t), 



TS(t) = £ 



> A ^ 2 



U e -2nTkt 



A > 
fc=l \l-KT I 



k 2 



2ttT 



(F9) 
2- (F10) 



For numerical evaluation purposes, it is useful to express 
this thermal sum in terms of Lerch transcendent func- 
tions: 



TS(t) = Re 



■#1 



A 



7o+«^ , 
2ttT ' 



27rTi 



s Ya 



<5 



1 \ 2ttT 
2 



2TrTt) 



(Fll) 



with the definitions of <?i(z;A), which is related to the 
Lerch transcendent function by $%{z; A) = $(e~ A , 1, z) — 
1/z, and of the symmetric part being 



Sy 2 [/(«) 



^Afe 



fe=l 

/(+z) + /(-z) 



(F12) 
(F13) 



The decay function is such that at the initial time it 
causes cancelation with the asymptotic values and the 
diffusion coefficients vanish. In this (asymptotic) high 
temperature perspective, the decay function contains two 
terms. The first decays at a cut-off dependent rate and 
can be expressed in closed form. The second decays with 
primarily temperature dependent rates and cannot be 
expressed in closed form with intuitive functions. It con- 
tains the initial time cancelation of the log(A/Q) diver- 
gence. Although well convergent at moderate times, the 
sum's contribution to the regular diffusion coefficient is 
very slow to converge at the initial time, even for mod- 
erate temperatures; see Fig. 14. 

While our expressions (F3)-(F4) can easily give us the 
zero temperature diffusion coefficients at asymptotically 
late time, they cannot easily give us the correspond- 
ing moderate time behavior in closed form. Moreover, 
the zero temperature limit of coth(w/2T) — > sgn(cj) 
means that our diffusion coefficient integrals cannot be 
cast as closed contour integrals. Nevertheless, the fre- 
quency integrals can be performed and the results ex- 
pressed in terms of exponential integrals with predictable 
time scales. At zero temperature (and in the high cut-off 
limit) we find the decay function to take the following 
form: 



lim DF(t) 

T->0 



(F14) 








<) 







Ei (At) 



A< 



-M 



where Ei(z) is the exponential integral, defined in A 2, 
which behaves like e~ z /z for large z. It should be noted 
that unlike the asymptotic limits of the diffusion coeffi- 
cients, the full time behavior is highly sensitive to the 
form of the cut-off regulator at low temperature. For our 
smooth regulator, we find relatively smoothly evolving 
diffusion coefficients (similar to the result in Ref. [17] at 
T = 1 Q) all the way down to zero temperature. In con- 
trast, a sharp cut-off of the form cx 9(u> — A) would 
produce the same average behavior, but with a slowly de- 
caying envelope modulating of considerable oscillations 
at the cut-off frequency. 

Analogous functions appear when we approximate the 
thermal sum in (FIO) [together with the first term on 
the right-hand side of (F9), which cancels any spurious 
poles at A = 2nTk] as an integral with a comparably soft 
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cut-off: 
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where fcj « 1. Still in the high cut-off limit, we find this 
qualitative approximation of the decay function to be 



DF(t) 



2 d 



Re 



7T dt 

2 d 

7T dt 



Ei( 27rTfej+7o+«0 t 



«o e -( 7o+ ^) t 

Ei([27rTfc 4 +zA] t) 



-iAt 



(F16) 



where we have discarded all finite terms at the initial 
time which decay at cut-off rates, as our approxima- 
tion ultimately ruins the behavior of DF(t) there. Thus, 
when using this approximate decay function, the time- 
dependent, decaying part of the diffusion coefficients 
must be "clamped" at the initial time. At moderate 
times, our approximation reveals the exact same form 
of exponential integral behavior as in the zero tempera- 
ture limit. But the temperature enters in such a way that 
the exponential decay inherent in Ei is not balanced out 
with a e _27rTfcit factor. Therefore, temperature is an in- 
herently stronger relaxation scale here [although there 
are additional e -70 ' factors from Gr(£) functions in the 
full diffusion coefficients]. 
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